Cholesky decomposition for S

[certik]

Does not work yet:

 SCF iteration:           1
 dsygvd returned info =         360
 The leading minor of order          120 of B is not positive definite. The factorization of B could not be completed and no eigenvalues or eigenvectors were computed.
Note: The following floating-point exceptions are signalling: IEEE_INVALID_FLAG
ERROR STOP eigh: dsygvd error

[HaoZeke]

Version: 1.2.2
Source dir: /home/rgoswami/Git/Github/Fortran/featom
Build dir: /home/rgoswami/Git/Github/Fortran/featom/bbdir
Build type: native build
Project name: featom
Project version: 0.1.0
Fortran compiler for the host machine: /home/rgoswami/micromamba/envs/fe/bin/ifort (intel 2021.6.0 "ifort (IFORT) 2021.6.0 20220226")
Fortran linker for the host machine: /home/rgoswami/micromamba/envs/fe/bin/ifort ld.bfd 2.41.0
Host machine cpu family: x86_64
Host machine cpu: x86_64
Found pkg-config: /home/rgoswami/micromamba/envs/fe/bin/pkg-config (0.29.2)
Run-time dependency mkl-dynamic-lp64-seq found: YES 2023.2
Build targets in project: 10

Generalized eigh:

With this PR, chol:

[HaoZeke]

! Replacing   D(:,:7) = matmul(invS(:,:,kappa), D2(:,:7))
        call dgemm('N', 'N', size(D,1), 7, size(D,1), 1.0d0, invS(:,:,kappa), size(D,1), D2(:,:7), size(D2,1), 0.0d0, D(:,:7), size(D,1))

[certik]

One should use:

         IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN
*
*           For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
*           backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y
*
            IF( UPPER ) THEN
               TRANS = 'N'
            ELSE
               TRANS = 'T'
            END IF
*
            CALL DTRSM( 'Left', UPLO, TRANS, 'Non-unit', N, M, ONE, B,
     $                  LDB, Z, LDZ )
*

[HaoZeke]

With DTRSM:

[certik]

--- a/tests/rlda/conv_rlda.f90
+++ b/tests/rlda/conv_rlda.f90
@@ -32,11 +32,11 @@ real(dp), allocatable :: R(:), Rp(:), V_tot(:), density(:)
 real(dp) :: eps
 integer :: p
 
-do p = 3, 8
+do p = 6, 6
     eps = 10.0_dp**(-p)
     eps = eps * 1.2_dp ! Allow numerical differences across compilers/platforms
     print *, "Test eps:", eps
-    do Z = 92, 1, -1
+    do Z = 92, 92, -1
         call get_LDA_energies(Z, ks_energies_exact, E_tot_exact)
         n_orb = size(ks_energies_exact)
         NN = get_N(Z, p)

Build and run using:

make -f Makefile.manual
cd tests/rlda
time ./conv_rlda

[HaoZeke]

dftatom comparison:
gfortran:

intel:

[HaoZeke]

With:

FC=$(which ifort) FFLAGS="-O3 -xHost -ipo -no-prec-div -fp-model fast=2 -profile-functions -g" meson setup bbdir -Dwith_tests=True --buildtype="release"
 vtune -collect hotspots ./bbdir/testDftDiracFast
    CPU Time: 0.440s
        Effective Time: 0.440s
        Spin Time: 0s
        Overhead Time: 0s
    Total Thread Count: 1
    Paused Time: 0s

Top Hotspots
Function                  Module            CPU Time  % of CPU Time(%)
------------------------  ----------------  --------  ----------------
[MKL LAPACK]@dsyevx       libmkl_core.so.2    0.160s             36.4%
assemble_radial_dirac_sh  libfeatom.so        0.100s             22.7%
[MKL LAPACK]@dsygst       libmkl_core.so.2    0.060s             13.6%
phih                      libfeatom.so        0.060s             13.6%
MKL_Load_Lib_Ex           libmkl_core.so.2    0.020s              4.5%
[Others]                  N/A                 0.040s              9.1%

[certik]

TODO:

• assembly (don't do S, use symmetry, no checks at the end)
• phih / dphi --- reduce usage, precalculate
• f2quad -- precalculate and do not call for GJ