ground state calculation with Lanczos and OpenMP

[jurekokalj]

Hi everyone,

as encouraged by Alex I would like to report a misbehaviour of ground state calculation with Lanczos (e.g. with eig0(...) ) that appears when using OpenMP and larger number of Lanczos steps.

The ground state is calculated in two steps to avoid high memory usage. In first step one runs Lanczos procedure, builds Lanczos (Krylov) vectors, builds tridiagonal matrix, which is then diagonalised. This gives the ground state energy and the ground state in terms of the Lanczos states's coefficients.
These coefficient are used in the next step in which one runs the same Lanczos procedure from the same starting vector. The ground state is gradually build from the corresponding coefficient and the Lanczos vector currently in memory. So, storing of all Lanczos vectors is not needed.

For this to work, Lanczos states must be exactly the same in both Lanczos runs and such is the case in single threaded no-openmp calculation and even for the ghost states. However, when openmp is used, the vectors are not exactly the same due to different numerical error buildup. This originates in different order of mathematical operations in Lanczos vector calculations due to many threads.

The calculated ground state therefore starts to diverge from the true ground state and the divergence is bigger for more lanczos steps, more openmp threads, more dynamical (vs. statical) openmp scheduling.

I believe this is hard to fix, since it would require to somehow force the same order of operations in both Lanczos procedures within multithreaded calculations. However, any ideas are welcome and maybe this report can be useful on its own.

For now it is good to stick with a lower number of Lanczos steps and maybe to testing the obtained ground state.

Best
Jure

[awietek]

Hi Jure,

thank you for reporting this and welcome to the XDiag discussion forum! As we discussed, I do think XDiag version 0.5 now mostly does static scheduling, but if you link to a Linear algebra backend, say an Intel MKL (which is called for dot products), then that package might have dynamic scheduling in it's parallelization. Could you give a minimal example together with which Linear Algebra library you used which show the behavior? One solution indeed could be to implement a dot product, which is statically scheduled in XDiag exactly for this purpose and I'd be willing to look into the issue.

Best,
Alex

[jurekokalj]

Hi Alex, Thanks for the reply. I attach a t-J example, which produced the output: J: 3.00000000e-01 T: 1.00000000e+00 block: tJ: nsites : 16 nup : 8 ndn : 6 irrep : defined with ID accebb8c026878c2 dimension: 22,470 ID : 230794386f1fc628 Number of omp threads 1 L.steps: 500. E0 form (Lanczos) and from <gs|H|gs>: -9.136796977749373 -9.136796977749377 L.steps: 500. E0 form (Lanczos) and from <gs|H|gs>: -9.136796977749373 -9.136796977749377 Number of omp threads 10 L.steps: 500. E0 form (Lanczos) and from <gs|H|gs>: -9.136796977749382 -9.13679697774668 L.steps: 500. E0 form (Lanczos) and from <gs|H|gs>: -9.136796977749386 -9.054268816594773 L.steps: 500. E0 form (Lanczos) and from <gs|H|gs>: -9.136796977749373 -9.04236454758743 L.steps: 500. E0 form (Lanczos) and from <gs|H|gs>: -9.136796977749396 -9.136796977739793 real   1m42.477s user   2m58.651s sys    1m2.298s This was calculated on my laptop with OSX, so there is no MKL used. BLAS is used. Additionally I set : export OMP_NUM_THREADS=10 export MKL_NUM_THREADS=1 export OPENBLAS_NUM_THREADS=1 so I believe the origin of misbehavior is not in MKL, but in OMP as well. Static scheduling does help as the FOR loops are equally distributed to OMP threads for both steps or both Lanczos runs. But different numerical error buildup also comes when atomic operation of calculated vector takes place, and this depends whether e.g. thread 1 updates the vector element first and then thread 2, or vice versa. This seems hard to control or to make it identical for both Lanczos runs. But this would be needed to have exactly the same Lanczos vectors in both steps. With many lanczos steps one gets ghost states - extra ground states. And the final accuracy of the ground states now depends on weather it has more weight in initial lanczos states (more similar for both runs) or later lanczos states (more different). Hope this is of any help. Best Jure From: Alexander Wietek ***@***.***> Reply to: awietek/xdiag ***@***.***> Date: Wednesday, 16 September 2026 at 10:00 To: awietek/xdiag ***@***.***> Cc: jurekokalj ***@***.***>, Author ***@***.***> Subject: Re: [awietek/xdiag] ground state calculation with Lanczos and OpenMP (Discussion #119) Hi Jure, thank you for reporting this and welcome to the XDiag discussion forum! As we discussed, I do think XDiag version 0.5 now mostly does static scheduling, but if you link to a Linear algebra backend, say an Intel MKL (which is called for dot products), then that package might have dynamic scheduling in it's parallelization. Could you give a minimal example together with which Linear Algebra library you used which show the behavior? One solution indeed could be to implement a dot product, which is statically scheduled in XDiag exactly for this purpose and I'd be willing to look into the issue. Best, Alex — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications, keep track of coding agent tasks and review pull requests on the go with GitHub Mobile for iOS and Android. Download it today! You are receiving this because you authored the thread.Message ID: ***@***.***>

[awietek]

I see, so you probably artificially increase the number of Lanczos iterations (500) beyond the convergence threshold to do e.g. FTLM, is that correct? And then after 500 steps things start to diverge. Okay, that's a good thing to check out, I will look into it.

[jurekokalj]

Yes. In some instances one would prefer number of Lanczos steps beyond convergence of ground state. For example, higher frequency resolution within FTLM, or better converge eigenstate in the middle of spectra for microcanonical Lanczos. From: Alexander Wietek ***@***.***> Reply to: awietek/xdiag ***@***.***> Date: Thursday, 17 September 2026 at 15:02 To: awietek/xdiag ***@***.***> Cc: jurekokalj ***@***.***>, Author ***@***.***> Subject: Re: [awietek/xdiag] ground state calculation with Lanczos and OpenMP (Discussion #119) I see, so you probably artificially increase the number of Lanczos iterations (500) beyond the convergence threshold to do e.g. FTLM, is that correct? And then after 500 steps things start to diverge. Okay, that's a good thing to check out, I will look into it. — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications, keep track of coding agent tasks and review pull requests on the go with GitHub Mobile for iOS and Android. Download it today! You are receiving this because you authored the thread.Message ID: ***@***.***>

[awietek]

Hi Jure,

I now looked into it and a think I solved the issue. You can check out the new branch "deterministic-apply", recompile the XDiag library, recompile your example and see whether this solves your issue.

The apply function is being used in a deterministic way, to a degree where even different thread counts should give you a bitwise identical result. What is nice about this, is that in the process I did get rid of some atomic operations, which now turns out to speed up all Lanczos runs by roughly 10%, so an extremely pleasant side effect.

Let me know if this works for you or if you need help with checking out the "deterministic-apply" branch. I would love to hear your feedback and whether it works for you. If so, I think this warrants a new minor release, since now basically all Lanczos runs are deterministic and have a 10% speedup, which is great.

Best,

Alex

[jurekokalj]

Hi Alex, This is amazing. It works. How did you do it and how did you make it deterministic for OpenMP runs – namely with the same order of operations. Is there some payoff, like higher memory usage? I did a short test. It took me a while, due to other workloads, but also since the runtime was a bit strange and ran it on laptop and compute node. Sorry for the delay. I ran 18 sites 2D t-J model with wavevector =(0,0). Attached it the code and xmol file, just in case. Maybe the system is still a bit small... The OMP part is with 16 threads. │ nsites   : 18 │ d        : 3 │ nup      : 9 │ ndn      : 7 │ permutation symmetries used │ irrep ID : b1463ef130171114 │ basis    : BasistJSymmetric<LinTable<uint32_t> > │ ID       : 8bab7d5e1508dc88 │ dimension: 97.240 Frist some old version of xdiag with ground state problem: XDiag Version: 0.4.1 Number of omp threads 1 L.steps: 400. E0 from (Lanczos) and from  <gs|H|gs>: -10.106700812385052 -10.106728044763393 30.25095 secs L.steps: 400. E0 from (Lanczos) and from  <gs|H|gs>: -10.106700812385052 -10.106728044763393 29.87334 secs Number of omp threads 16 L.steps: 400. E0 from (Lanczos) and from  <gs|H|gs>: -10.106700812385041 -9.861265985267721 9.71205 secs L.steps: 400. E0 from (Lanczos) and from  <gs|H|gs>: -10.106700812385048 -10.101890008093186 9.60589 secs L.steps: 400. E0 from (Lanczos) and from  <gs|H|gs>: -10.106700812385052 -10.007881015500523 9.60014 secs L.steps: 400. E0 from (Lanczos) and from  <gs|H|gs>: -10.106700812385043 -9.742908525385396 9.60702 secs Then the deterministic version, which nicely gives the correct gs and is a bit slower than the old version. XDiag Version: 0.5.1 (deterministic) Number of omp threads 1 L.steps: 400. E0 from (Lanczos) and from  <gs|H|gs>: -10.10670081238505 -10.106603102046439 75.80924 secs L.steps: 400. E0 from (Lanczos) and from  <gs|H|gs>: -10.10670081238505 -10.106603102046439 75.80656 secs Number of omp threads 16 L.steps: 400. E0 from (Lanczos) and from  <gs|H|gs>: -10.10670081238505 -10.106603102046439 17.49550 secs L.steps: 400. E0 from (Lanczos) and from  <gs|H|gs>: -10.10670081238505 -10.106603102046439 17.50237 secs L.steps: 400. E0 from (Lanczos) and from  <gs|H|gs>: -10.10670081238505 -10.106603102046439 17.50573 secs L.steps: 400. E0 from (Lanczos) and from  <gs|H|gs>: -10.10670081238505 -10.106603102046439 17.45836 secs And last the newest github version.  Here I am not sure, why it is so slow. Seem to be slow on my laptop OSX and on compute nodes (I guess rocky linux 8.8). Possibly I am doing something wrong. XDiag Version: 0.5.1 Number of omp threads 1 L.steps: 400. E0 from (Lanczos) and from  <gs|H|gs>: -10.10670081238505 -10.106686905600789 699.11764 secs L.steps: 400. E0 from (Lanczos) and from  <gs|H|gs>: -10.10670081238505 -10.106686905600789 698.81036 secs Number of omp threads 16 L.steps: 400. E0 from (Lanczos) and from  <gs|H|gs>: -10.106700812385032 -9.838697961000893 644.11569 secs L.steps: 400. E0 from (Lanczos) and from  <gs|H|gs>: -10.106700812385045 -10.10632728537254 642.40128 secs L.steps: 400. E0 from (Lanczos) and from  <gs|H|gs>: -10.106700812385046 -9.92478349220631 641.91537 secs L.steps: 400. E0 from (Lanczos) and from  <gs|H|gs>: -10.106700812385037 -9.668513641342523 641.55887 secs Anyhow. The deterministic version seems to work great. It is with my compilation a bit slower then the older version, but is definitely very useful. Thanks. Best Jure From: Alexander Wietek ***@***.***> Reply to: awietek/xdiag ***@***.***> Date: Wednesday, 23 September 2026 at 10:20 To: awietek/xdiag ***@***.***> Cc: jurekokalj ***@***.***>, Author ***@***.***> Subject: Re: [awietek/xdiag] ground state calculation with Lanczos and OpenMP (Discussion #119) Hi Jure, I now looked into it and a think I solved the issue. You can check out the new branch "deterministic-apply", recompile the XDiag library, recompile your example and see whether this solves your issue. The apply function is being used in a deterministic way, to a degree where even different thread counts should give you a bitwise identical result. What is nice about this, is that in the process I did get rid of some atomic operations, which now turns out to speed up all Lanczos runs by roughly 10%, so an extremely pleasant side effect. Let me know if this works for you or if you need help with checking out the "deterministic-apply" branch. I would love to hear your feedback and whether it works for you. If so, I think this warrants a new minor release, since now basically all Lanczos runs are deterministic and have a 10% speedup, which is great. Best, Alex — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications, keep track of coding agent tasks and review pull requests on the go with GitHub Mobile for iOS and Android. Download it today! You are receiving this because you authored the thread.Message ID: ***@***.***>

[awietek]

Uh, this looks like quite a bad performance issue. Could you share the exact Hamiltonian and sector you are running this with?

[jurekokalj]

Sure. Please find attached .cpp and .toml.   Best Jure From: Alexander Wietek ***@***.***> Reply to: awietek/xdiag ***@***.***> Date: Wednesday, 7 October 2026 at 20:44 To: awietek/xdiag ***@***.***> Cc: jurekokalj ***@***.***>, Author ***@***.***> Subject: Re: [awietek/xdiag] ground state calculation with Lanczos and OpenMP (Discussion #119) Uh, this looks like quite a bad performance issue. Could you share the exact Hamiltonian and sector you are running this with? — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications, keep track of coding agent tasks and review pull requests on the go with GitHub Mobile for iOS and Android. Download it today! You are receiving this because you authored the thread.Message ID: ***@***.***>

[awietek]

Hi Jure, the attachment seems to be missing, can you just send them to me via email?

[jurekokalj]

Sorry. Maybe replying via email removes the attachments. This should work. From: Alexander Wietek ***@***.***> Reply to: awietek/xdiag ***@***.***> Date: Thursday, 8 October 2026 at 09:08 To: awietek/xdiag ***@***.***> Cc: jurekokalj ***@***.***>, Author ***@***.***> Subject: Re: [awietek/xdiag] ground state calculation with Lanczos and OpenMP (Discussion #119) Hi Jure, the attachment seems to be missing, can you just send them to me via email? — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications, keep track of coding agent tasks and review pull requests on the go with GitHub Mobile for iOS and Android. Download it today! You are receiving this because you authored the thread.Message ID: ***@***.***>

[awietek]

Hi Jure,

I just ran the benchmarks on the latest main branch and did not find a performance regression, the code is actually faster than with 0.4.1. There was an intermediate version, that might explain your observed behavior. Hence, I would just like to ask you to check out the latest code on the main branch and run your calculation again. Please let me know if that works for you.

All the best,
Alex

[jurekokalj]

Hi Alex, Yes that worked. I was compiling with older gcc and I guess that was the reason for regression. Now with the newest version 0.5.1 and by using newer modules  foss/2026.1, the calculation is much faster. Here is the output of the test. │ nsites : 18 │ d : 3 │ nup : 9 │ ndn : 7 │ permutation symmetries used │ irrep ID : b1463ef130171114 │ basis : BasistJSymmetric<LinTable<uint32_t> > │ ID : 8bab7d5e1508dc88 │ dimension: 97.240 Number of omp threads 1 L.steps: 400. E0 from (Lanczos) and from <gs|H|gs>: -10.106700812385043 -10.107671090341684 65.83691 secs L.steps: 400. E0 from (Lanczos) and from <gs|H|gs>: -10.106700812385043 -10.107671090341684 65.75003 secs Number of omp threads 16 L.steps: 400. E0 from (Lanczos) and from <gs|H|gs>: -10.106700812385043 -10.107671090341698 16.33685 secs L.steps: 400. E0 from (Lanczos) and from <gs|H|gs>: -10.106700812385043 -10.107671090341698 16.41696 secs L.steps: 400. E0 from (Lanczos) and from <gs|H|gs>: -10.106700812385043 -10.107671090341698 16.38487 secs L.steps: 400. E0 from (Lanczos) and from <gs|H|gs>: -10.106700812385043 -10.107671090341698 16.41520 secs It works very fine. Thanks. Jure From: Alexander Wietek ***@***.***> Reply to: awietek/xdiag ***@***.***> Date: Friday, 9 October 2026 at 13:01 To: awietek/xdiag ***@***.***> Cc: jurekokalj ***@***.***>, Author ***@***.***> Subject: Re: [awietek/xdiag] ground state calculation with Lanczos and OpenMP (Discussion #119) Hi Jure, I just ran the benchmarks on the latest main branch and did not find a performance regression, the code is actually faster than with 0.4.1. There was an intermediate version, that might explain your observed behavior. Hence, I would just like to ask you to check out the latest code on the main branch and run your calculation again. Please let me know if that works for you. All the best, Alex — Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications, keep track of coding agent tasks and review pull requests on the go with GitHub Mobile for iOS and Android. Download it today! You are receiving this because you authored the thread.Message ID: ***@***.***>