Slepc with preconditionig converges very fast to the wrong result

sangallidavide commented 5 months ago

The bug can be easely reproduced using the following input for C6H3Cl3 and comparing with 02_ALDA_diago_ws

N.B.: you need to use the "SAVE_converted" folder with recent versions of yambo

mv SAVE SAVE_old
mv SAVE_converted SAVE
yambo
yambo -F bse.in

with the following bse.in

optics                       # [R OPT] Optics
bse                          # [R BSE] Bethe Salpeter Equation.
tddft                        # [R   K] Use TDDFT kernel
bsk                          # [R BSK] Bethe Salpeter Equation kernel
bss                          # [R BSS] Bethe Salpeter Equation solver
NonPDirs= "XYZ"              # [X/BSS] Non periodic chartesian directions (X,Y,Z,XY...)
BSEmod= "resonant"           # [BSE] resonant/causal/coupling
BSKmod= "ALDA"               # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
BSSmod= "s"                  # [BSS] (h)aydock/(d)iagonalization/(i)nversion/(t)ddft`
BSENGexx= 1535         mHa   # [BSK] Exchange components

BSSNEig= 30                     # [SLEPC] Number of eigenvalues to compute
BSSEnTarget= 0.000000      eV    # [SLEPC] Target energy to find eigenvalues
BSSSlepcApproach= "Krylov-Schur" # [SLEPC] Approach ("Krylov-Schur","Generalized-Davidson","Jacob-Davidson")
BSSSlepcPrecondition= "preonly+jacobi" # [SLEPC] Precondition technique (none|preonly+jacobi|bcgs+jacobi)
BSSSlepcExtraction= "ritz"       # [SLEPC] Extraction technique (ritz|harmonic)
BSSSlepcMaxIt=0                  # [SLEPC] Maximum number of iterations
BSSSlepcNCV=0                    # [SLEPC] Dimension of the subspace
BSSSlepcTol= 0.100000E-5         # [SLEPC] Tolerance for the iterative solver
BSSSlepcMatrix                # [SLEPC] Store slepc matrix, faster but more memory consuming

rim_cut
CUTGeo= "ws xyz"             # [CUT] Coulomb Cutoff geometry: box/cylinder/sphere
CUTwsGvec= 0.700000            # [CUT] WS cutoff: number of G to be modified

% BEnRange
  3.00000 | 11.00000 | eV    # [BSS] Energy range
%
% BDmRange
  0.10000 |  0.10000 | eV    # [BSS] Damping range
%
BEnSteps=  800               # [BSS] Energy steps
% BLongDir
 1.000000 | 1.000000 | 1.000000 |        # [BSS] [cc] Electric Field
%
% BSEBands
  14 |  34 |                 # [BSK] Bands range
%

Also in the code I tried to play with different options. Relevant info in file src/linear_algebra/MATRIX_slepc.F

 !
 ! * Krylov subspaces: EPSKRYLOVSCHUR, STPPRECOND not accepted,
 !                        STSINVERT + KSPPREONLY+ PCJACOBI gives wrong eigenvalues. Too strong (?)
 !                                  +  KSPBCGS  + PCJACOBI works fine
 !
 ! * Generalized-Davidson: EPSGD, RIGHT VECTORS ONLY
 !                       STPPRECOND + KSPPREONLY + PCJACOBY very fast for few eigenvectors
 !
 ! * Jacob-Davidson: EPSJD, NOT WORKING WITH KSP PREONLY
 !                    STPPRECOND + KSPCG SUPER-SLOW + PCJACOBY super slow
 !

But I was not sure which was the correct set of combinations to use

sangallidavide commented 5 months ago

Results with (wrong) and without (ok) preconditioning

log with preconditioning

<04s> [SLEPC] Approach                          : Krylov-Schur
<04s> [SLEPC] Precondition method               : preonly+jacobi
<04s> [SLEPC] Extraction method                 : ritz
<04s> [SLEPC] Number of requested eigenvalues   :  30
<04s> [SLEPC] Criterion is target energy        :  0.000000 [eV]
<04s> [SLEPC] Stopping condition tolerance      :  0.100000E-5
<04s> [SLEPC] Stopping condition max iterations : -2
<04s> [SLEPC] Iteration #1 - converged States 5 - error :  0.185728E-8  0.192937E-5
<04s> [SLEPC] Iteration #2 - converged States 24 - error :  0.142434E-6  0.425348E-5
<04s> [SLEPC] Iteration #3 - converged States 38 - error :  0.373440E-8  0.443773E-4
<04s> [SLEPC] Number of iterations              :  3
<04s> [SLEPC] Number of eigenvalues        [NEV]:  30
<04s> [SLEPC] Max. subspace size of solver [NCV]:   60
<04s> [SLEPC] Max. allowed dim             [MPD]:   60
<04s> [SLEPC] Number of converged states        :   38
<04s> EPS R residuals |########################################| [100%] --(E) --(X)
<04s> BSK resp. funct |########################################| [100%] --(E) --(X)

log without preconditionig:

<04s> [SLEPC] Approach                          : Krylov-Schur
<04s> [SLEPC] Extraction method                 : ritz
<04s> [SLEPC] Number of requested eigenvalues   :  30
<04s> [SLEPC] Criterion is target energy        :  0.000000 [eV]
<04s> [SLEPC] Stopping condition tolerance      :  0.100000E-5
<04s> [SLEPC] Stopping condition max iterations : -2
<04s> [SLEPC] Iteration #1 - converged States 5 - error :  0.197869E-8  0.243883E-5
<04s> [SLEPC] Iteration #2 - converged States 9 - error :  0.453823E-6  0.125730E-5
<04s> [SLEPC] Iteration #3 - converged States 26 - error :  0.483774E-6  0.263769E-5
<04s> [SLEPC] Iteration #4 - converged States 38 - error :  0.103942E-6  0.334431E-4
<04s> [SLEPC] Number of iterations              :   4
<04s> [SLEPC] Number of eigenvalues        [NEV]:  30
<04s> [SLEPC] Max. subspace size of solver [NCV]:   60
<04s> [SLEPC] Max. allowed dim             [MPD]:   60
<04s> [SLEPC] Number of converged states        :   38
<04s> EPS R residuals |########################################| [100%] --(E) --(X)
<04s> BSK resp. funct |########################################| [100%] --(E) --(X)

joseeroman commented 5 months ago

The combination BSSSlepcApproach="Krylov-Schur" with BSSSlepcPrecondition="preonly+jacobi" should be forbidden.

Generalized Davidson and Jacobi-Davidson (by the way, it is misspelled as Jacob-Davidson throughout yambo) are preconditioned eigensolvers, which means that they can optionally apply a preconditioner (PC in PETSc) at each iteration of the eigensolver. A good preconditioner may improve convergence substantially. In the case of Jacobi-Davidson, the solver uses a combination KSP+PC, where KSP is an iterative linear solver, so instead of just applying a preconditioner it "solves" a system with low accuracy (in practice, with a limited number of iterations). In this case one could set KSP=preonly to indicate "preconditioner only", but that would prevent the method from progressing in most cases, so it is not recommended.

In contrast, Krylov-Schur is a Krylov method, i.e., not a preconditioned method. This implies that when using shift-and-invert to compute eigenvalues close to a target, the associated linear systems must be solved accurately, with a full KSP+PC, and hence KSP must be something like bcgs or gmres, but not preonly.

Regarding Generalized Davidson and Jacobi-Davidson:

Generalized Davidson with diagonal preconditioner (jacobi) will be effective if the matrix is strongly diagonally dominant, which is probably the case in yambo.
Jacobi-Davidson often does not improve on Generalized Davidson if the latter works well. Also, the use of KSP=cg may not be appropriate in Jacobi-Davidson because CG requires a symmetric (Hermitian) matrix but Jacobi-Davidson modifies the input matrix in a way that symmetry may be lost.

sangallidavide commented 5 months ago

Ok. So I understand.

1) We should put a barrier for preconditioning with Krylov-Schur

2) We may consider using more Generalized-Davidson instead of always sticking to Krylov-Schur (although only in some cases)

Is that right?

sangallidavide commented 4 months ago

I'm checking now the connected code. But I'm not sure on how to change the following:

 if (l_precondition) then                                               
   !                                                                    
   if (epskind==EPSKRYLOVSCHUR) stkind=STSINVERT                        
   if (epskind==EPSGD)          stkind=STPRECOND                        
   if (epskind==EPSJD)          stkind=STPRECOND                        
   !                                                                    
   ! Default                                                            
   if (stkind==STSINVERT)                           kspkind=KSPBCGS     
   if (stkind==STPRECOND .and. epskind==EPSGD)      kspkind=KSPPREONLY  
   if (stkind==STPRECOND .and. epskind==EPSJD)      kspkind=KSPBCGS     
   ! From user                                                          
   if(index(BSS_slepc_precondition,KSPPREONLY)/=0)  kspkind=KSPPREONLY  
   if(index(BSS_slepc_precondition,KSPBCGS)/=0)     kspkind=KSPBCGS     
   !                                                                    
   ! Default                                                            
   pckind=PCJACOBI                                                      
   ! From user                                                          
   if(index(BSS_slepc_precondition,PCJACOBI)/=0)    pckind=PCJACOBI

Suggestions are welcome

andreamarini commented 4 months ago

Dear all. I would suggest to move this discussion to the https://github.com/yambo-code/yambo-libraries repo.

This issue is indeed related to a specific library.

yambo-code / yambo

Slepc with preconditionig converges very fast to the wrong result #68