Sorry, I'll disagree. When you write high performance code like this, there's no artistic expression, but different approaches.
These algorithms' design starts on paper, as pure maths. You get the analytic form, and convert this to numerical/discrete form and start to design from there.
You need to add normalization, subnormal number prevention, maybe some architecture specific optimization if you are targeting a specific CPU generation, etc. Then it gets optimized further as you run tests, profile the code and run your test cases.
So, no BLAS package looks the same, because they make different trade-offs, have different programming models, and handle things differently, at different points. Even more so, a newer generation processor with much higher FPU performance can have lower performance than previous generation, so you implement another copy of that function or add to detect that processor and have different paths just for that one.
IOW, when doing numerical maths on processors, esp. if you're doing this on many machines with MPI or multiple cores with hand rolled code or OpenMP, there's no single way. Maybe a best way for the specific CPU you have, but it still depends on the trade-offs and design choices you made in the very beginning.
These algorithms' design starts on paper, as pure maths. You get the analytic form, and convert this to numerical/discrete form and start to design from there.
You need to add normalization, subnormal number prevention, maybe some architecture specific optimization if you are targeting a specific CPU generation, etc. Then it gets optimized further as you run tests, profile the code and run your test cases.
So, no BLAS package looks the same, because they make different trade-offs, have different programming models, and handle things differently, at different points. Even more so, a newer generation processor with much higher FPU performance can have lower performance than previous generation, so you implement another copy of that function or add to detect that processor and have different paths just for that one.
IOW, when doing numerical maths on processors, esp. if you're doing this on many machines with MPI or multiple cores with hand rolled code or OpenMP, there's no single way. Maybe a best way for the specific CPU you have, but it still depends on the trade-offs and design choices you made in the very beginning.