Single precision computations to improve speed

[philipturner]

Have the authors of this repo considered using FP32 for the majority of the computations, instead of FP64? FP64 is the de facto standard in scientific computing, but I've always been able to get by with FP32, provided proper considerations are taken. There are often opportunities to compute bulk integrals in FP32, then accumulate into FP64, or only perform critical parts of a program in FP64.

I believe this is appropriate and may be taken seriously, after finding multiple commits in the repo's history, about 2% or 30% gains in various places, due to small changes in the code. Mixed precision is a very important part of the design space that deserves consideration. I estimate a 1.5x improvement in whole-program execution time on vector ALUs where FP32 throughput is 2x that of FP64 (I don't know how Fortran vectorizes the code; I typically work in languages with explicit vector types like SIMD8<Float32>).

As context, I did an analysis of the bottlenecks contributing to latency for 1 singlepoint with GFN2-xTB. I found that linear algebra kernels from external libraries consume 70% of the time (provided the integrals are properly parallelized with OpenMP). I am attempting to link xTB to a custom linear algebra library that casts the FP64 inputs to FP32 and does eigendecomposition in FP32, with faster speed. I expect issues, but with careful analysis the bugs or threshold violations can be fixed, making the program work.

https://github.com/philipturner/swift-xtb/blob/f887f2ca1dfad9c2c284764eb6377a7889004f60/Documentation/Archive/Experiments/OptimizingDiagonalization/OptimizingDiagonalization.md

However, I cannot rewrite the entire xTB codebase from scratch just to get the integrals ported to FP32. I expect that once dsyevd_ ssyevd_ is vastly sped up, integrals may be a non-negligible contributor to latency. In any case, every component of the program should be considered for latency reduction.

[foxtran]

about 2% or 30% gains in various places, due to small changes in the code.

There are extra 3-5% gains from not published yet PRs. And extra 10% gains can be achieved with runtime changes which are rarely using.

integrals may be a non-negligible contributor to latency.

Integrals are evaluating once per each geometry optimization step, so I do not think that it will help a lot.

every component of the program should be considered for latency reduction.

There are several issues regarding this. Once you decrease latency for quite big systems, the code becomes insufficient for small systems and that is a big problem since sometimes guys run millions of small molecules with xTB.

[foxtran]

I am attempting to link xTB to a custom linear algebra library that casts the FP64 inputs to FP32 and does eigendecomposition in FP32, with faster speed.

You can introduce FP32 solver here for first iterations while abs(eel - eold) > 1e-6. It should be quite safe.

xtb/src/scc_core.f90

Line 475 in d3f212b

call solver%fact_solve(env, H, S_factorized, emo)

[philipturner]

There are several issues regarding this. Once you decrease latency for quite big systems, the code becomes insufficient for small systems and that is a big problem since sometimes guys run millions of small molecules with xTB.

I'm targeting both small and large systems. Especially, my use case often involves long-time MD simulations with systems having 400 orbitals or less. Up to about ~1000 orbitals, latency scales with the square of orbital count, not the cube. It is here that integrals (quadratic scaling) have cost comparable to eigendecomposition (also quadratic scaling here).

You can introduce FP32 solver here for first iterations while abs(eel - eold) > 1e-6. It should be quite safe.

Thanks. I expect that, since the codebase is mostly in FP64, there are several error thresholds tripped by small changes sensitive to the shift from FP64 to FP32. Especially the "set numerical accuracy" (???) function from the C API.

[philipturner]

Integrals are evaluating once per each geometry optimization step, so I do not think that it will help a lot.

Perhaps it's something else contributing to latency. I've measured the time taken to compute dsyevd_, multiplied it by the SCF iterations per singlepoint, and the numbers didn't add up. However, I have Xcode profiler data clearly showing that at least 20% of the execution time per SCF cycle comes from stuff other than linear algebra calls.

Gradients are what's evaluated once per geometry optimization step, based on my understanding. For each SCF iteration, you construct the Hamiltonian with semiempirical approximations to the true real-space Poisson equation and similar. Then diagonalize the Hamiltonian and repeat until the residual error converges. Is this correct?

[philipturner]

Thanks. I expect that, since the codebase is mostly in FP64, there are several error thresholds tripped by small changes sensitive to the shift from FP64 to FP32. Especially the "set numerical accuracy" (???) function from the C API.

Are the majority or the minority of the SCF iterations in the region where the error exceeds 1e-6? Ideally, 100% of the iterations are in this region, and the calculation has fully converged once it got so low. Note that 1e-6 is still far above the dynamic range of FP32 (1e-38). The error itself can be summed in FP64 while the terms contributing to it can originate in FP32.

[foxtran]

Perhaps it's something else contributing to latency.

You can look named profile with timings with Tracy using #1284.

However, I have Xcode profiler data clearly showing that at least 20% of the execution time per SCF cycle comes from stuff other than linear algebra calls.

Ah... You are using Mac with gfortran builds of xtb... There is a big room for optimizing gfortran builds to match Intel compiler performance.

For each SCF iteration, you construct the Hamiltonian with semiempirical approximations to the true real-space Poisson equation and similar. Then diagonalize the Hamiltonian and repeat until the residual error converges. Is this correct?

Something like that. Some parts of Hamiltonian are cached and are computing once at the beginning of SCF iterations (see src/scf_module.f90).

Are the majority or the minority of the SCF iterations in the region where the error exceeds 1e-6? Ideally, 100% of the iterations are in this region, and the calculation has fully converged once it got so low. Note that 1e-6 is still far above the dynamic range of FP32 (1e-38). The error itself can be summed in FP64 while the terms contributing to it can originate in FP32.

No. The idea is to make first iterations cheaper while FP32 accuracy allows us to do it, then switch to FP64 precision.

[philipturner]

Something like that. Some parts of Hamiltonian are cached and are computing once at the beginning of SCF iterations (see src/scf_module.f90).

I don't have a lot of experience with how semiempirical TB works. Only, that I attempted to make a full ab initio simulator, based on the physical equations alone, with real-space adaptive mesh refinement and all-electron calculations. Got stuck on the kinetic energy operator not enjoying this basis. And it taking 1 minute (!) per singlepoint with linear-scaling DFT on supercomputers. Semiempirical looks much more complicated to code; I can't code the entire thing myself from the Hartree-Fock equation alone.

Ah... You are using Mac with gfortran builds of xtb... There is a big room for optimizing gfortran builds to match Intel compiler performance.

Mac (Apple silicon) has a dedicated linear algebra accelerator on their CPU, called the AMX (https://github.com/corsix/amx). It is especially well suited for linear algebra with small matrices. It improves FP64 matrix throughput by 2x over FP64 vector, but the beast is FP32 matrix throughput - 8x more than FP64 on ARM NEON vector units (4x more than FP32 vector). I've seen measurable gains from switching from OpenBLAS to Apple Accelerate, such as 3x latency reduction at ~190-360 orbitals. The benchmarks I logged were already using Accelerate, but I'm seeking to improve speed even further beyond that.

No. The idea is to make first iterations cheaper while FP32 accuracy allows us to do it, then switch to FP64 precision.

There's likely room for a "mixed precision" where FP32 is used for the less sensitive parts, and FP64 for the more sensitive. That should extend out the region where FP32 can produce benefits, more than just FP32 alone. Does that not seem reasonable?

---

There's a lot of misconceptions about FP32 being "inaccurate", which go uncorrected because the HPC community just solves it with supercluster hardware dedicated to FP64 primarily. Look at the offerings by AMD's CDNA architecture. Consumer hardware has focused more on FP32, both GPUs and novel matrix accelerators attached to the CPU. Apple in particular.

Even within the diagonalization kernel, there's regions that could be in FP32 while other critical regions stay in FP64. And I have the motivation to explore these techniques because Moore's Law is dead, and consumer hardware won't clock any faster.

Thanks for the insights - I might come off a little frustrated, but I care deeply about the topic.

[philipturner]

And it taking 1 minute (!) per singlepoint with linear-scaling DFT on supercomputers.

In that paper by Jean-Luc Fattebert, they just switched all of their data from FP64 to FP32 and only accumulated the sums of the ab-initio integrals in FP64. Which happens to be the bottleneck in ab initio simulations (matrix diagonalization is a minor fraction of execution latency).

Strangely, Fattebert had no problems. But it's never occurred to almost every real-space ab initio program I've surveyed, such as Octopus, to simply reassess the need for 100% of the calculations to be in FP64. Meanwhile, my attempts at an ab initio simulator deliberatedly did 100% of the calculations in FP32, and I had no problems.

[foxtran]

Mac (Apple silicon) has a dedicated linear algebra accelerator on their CPU

That is nice. What I'm saying that GCC produce worse x86_64 assembler in comparison with Intel compilers. So, binaries compiled with GCC are slower than binaries compiled with Intel compilers if they use the same math library.

Does that not seem reasonable?

It makes sense. One just needs to code it.

In that paper by Jean-Luc Fattebert, they just switched all of their data from FP64 to FP32 and only accumulated the sums of the ab-initio integrals in FP64.

Ab-Initio integrals should have 11-12 significant digits, so, it may work.

Meanwhile, my attempts at an ab initio simulator deliberatedly did 100% of the calculations in FP32, and I had no problems.

Here, in xtb, are several issues that some thresholds excludes some interactions (even if they are at 1e-11 level) that often breaks geometry optimizations. More specifically, you may find multiple geometries which are really close to each other (shift of single atom at 1e-5 A) but their energies will be different to up to several kcal/mol.

[philipturner]

This is all content that ought to be explored in more detail. I'm confident there's room for speed improvement, although not without a lot of effort.

That is nice. What I'm saying that GCC produce worse x86_64 assembler in comparison with Intel compilers. So, binaries compiled with GCC are slower than binaries compiled with Intel compilers if they use the same math library.

I understand. Apple and Intel are both very different hardware, and very different compiler ecosystems.

It makes sense. One just needs to code it.

I'm willing to take the time to code it.

Ab-Initio integrals should have 11-12 significant digits, so, it may work.

I strongly question the number 11-12 significant digits. Can you cite a peer-reviewed source?

Here, in xtb, are several issues that some thresholds excludes some interactions (even if they are at 1e-11 level) that often breaks geometry optimizations. More specifically, you may find multiple geometries which are really close to each other (shift of single atom at 1e-5 A) but their energies will be different to up to several kcal/mol.

My use case is interactive molecular dynamics simulations. Forces are much less sensitive to small things like changes in energy. Eric Drexler highlighted this in Nanosystems (1992) and Peter Eastman of OpenMM has a great paper:

Energy Conservation as a Measure of Simulation Accuracy

[philipturner]

I'm willing to take the time to code it.

There's actually a whole paper on this, where some Korean researchers used "dynamic precision" to speed up ab initio simulations ~8x on consumer RTX GPUs, where FP64 throughput is 1/32 of FP32. I can provide the source if you want.

[philipturner]

shift of single atom at 1e-5 A) but their energies will be different to up to several kcal/mol.

If a shift of 1 femtometer causes my structure to explode, I obviously have not engineered a very stable molecule for use in the real world.

[foxtran]

If a shift of 1 femtometer causes my structure to explode, I obviously have not engineered a very stable molecule for use in the real world.

It is not an explosion: it looks like Heaviside function.

[philipturner]

Perhaps a molecule prone to the Jahn-Teller distortion?

[philipturner]

You can look named profile with timings with Tracy using #1284.

As an aside, I was thinking of exactly the same thing. Had a Google Doc characterizing the latency/orbital^2 in microseconds/orbital^2. Lacked enough data to break down each linear algebra call's contribution.

I'll try this once I can get xTB compiling from source. I'm facing a bug where it crashes at systems with 770 orbitals or higher, meaning I can't even evaluate performance of the latest version against my own test cases. Learning Fortran from scratch at the moment, so I can start to reason about how in the world to compile this.