SLEPc solver with nest matrix format

[blmelp]

Joint work with @joseeroman

Introduces a new matrix format for the SLEPc solver, which optimizes the performance for the BSE with coupling. This format can be used only if BS_K_coupling and l_BS_ares_from_res are TRUE.

The new format uses a PETSc nest matrix to take advantage of the structure of the matrix in this case:

| R      C  |
|-C^*   -R^T|

The top blocks are stored explicitly, and the bottom blocks are shell matrices, reducing the memory required to store the matrix by approximately half.

Previous BSSSlepcMatrix logical option is changed to BSSSlepcMatrixFormat, with three options:

• 'shell': Shell matrix
• 'explicit': If it is posible, it will use the new format with a nest matrix to store the four blocks
• 'explicit-debug': It will always use the old version of the explicit matrix, without the nest storage

If Yambo and PETSc are configured with CUDA, matrix operations with the nest format will be performed in GPU.

[sangallidavide]

Just a general remark.

The standard structure of BSE is

       |   R      C    |                             
   K = |               |                  
       | -C^*   -R^*   |                             

This is equivalent to the above used expression since R is hermitian.

Notice that yambo uses a modified BSE matrix, with a modified coupling, D= i C It still holds the general structure

       |    R      D    |                             
   K = |                |                   
       | - D^*   -R^*   |                             

Please check that this does not have any impact on the implementation. The modified coupling is introduced to obtain a pseudo-hermitian matrix also in presence of fractional occupations. For more details please check eq 17 here: https://journals.aps.org/prb/abstract/10.1103/PhysRevB.93.195205

[blmelp]

I do not think so, because our change does not modify anything about the computation of the values of the coupling matrix, and the transformations from C to -C^ and from D to -D^ are the same, so it should not be an issue. @joseeroman what do you think?

[joseeroman]

The change from C to D=i C should be ok. Regarding the comment about R, I am confused because if R is hermitian then R^T is not the same as R^*, but R^*=conj(R^T). Am I missing something?

[sangallidavide]

The change from C to D=i C should be ok. Regarding the comment about R, I am confused because if R is hermitian then R^T is not the same as R^, but R^=conj(R^T). Am I missing something?

With R^* I mean conjg(R) I use R^+ for hermitian(R)

Latex style. Is this the source of confusion ?

[joseeroman]

Ok, everything clear now.

[sangallidavide]

I realized now that there is a bug which shows up when running the new nested approach in parallel.

Nested approach in serial works fine. Explicit-debug approach in parallel also works fine.

I can provide a test if needed. Just let me know.

[blmelp]

If you can, send us the test. We will look into it

[sangallidavide]

On lynx in the folder /home/sangalli/calcs/isolv_bse/bug-nested-approach you have everything you need to reproduce the bug

you can run the script reproducer.sh and compare the two output files vimdiff slepc1_p*/o-slepc1_p*.alpha*

[blmelp]

We have created pull request #94 to fix the bug but we are not sure how to test it. Can you give it a try?

[sangallidavide]

Checked, and it works. Pull request merged.

Can you remember me from which version of slepc/petsc the nest matrix format is available ?

[blmelp]

It is better to use 3.21 of both PETSc and SLEPc (they must be the same version).