Initial implementation of Cholesky

Author: certik

src/dirac.f90

@@ -16,6 +16,7 @@ module dirac
     nlsf2focc, get_atomic_states_rel, nlf2focc, get_atomic_states_nonrel
 use energies, only: thomas_fermi_potential
 use iso_c_binding, only: c_double, c_int
+use lapack, only: dpotrf
 implicit none
 private
 public solve_dirac, csolve_dirac, solve_dirac_eigenproblem
@@ -37,7 +38,8 @@ subroutine solve_dirac_eigenproblem(Nb, Nq, Lmin, Lmax, alpha, alpha_j, xe, xiq_
 real(dp), intent(out) :: V(:,:)
 integer, intent(out) :: idx
 real(dp), intent(in) :: E_dirac_shift
-integer :: kappa, i
+real(dp) :: SU(size(S,1),size(S,2))
+integer :: kappa, i, info, j
 idx = 0
 do kappa = Lmin, Lmax
     if (kappa == 0) cycle
@@ -58,6 +60,13 @@ subroutine solve_dirac_eigenproblem(Nb, Nq, Lmin, Lmax, alpha, alpha_j, xe, xiq_
     ! we still need two seperate eigensolves for 1e-8 accuracy:
     !H = (H + transpose(H))/2
     !S = (S + transpose(S))/2
+    SU = S
+    do j = 1, size(SU,1)
+        SU(j+1:,j) = 0
+    end do
+    call dpotrf('U', size(SU,1), SU, size(SU,1), info)
+    if (info /= 0) error stop
+    S = matmul(transpose(SU), SU)
 
     if (accurate_eigensolver) then
         call eigh(H, S, lam)

src/lapack.f90

@@ -89,6 +89,13 @@ SUBROUTINE DSYGVX( ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB, &
                        Z( LDZ, * )
     END SUBROUTINE
 
+    SUBROUTINE DPOTRF( UPLO, N, A, LDA, INFO )
+    import :: dp
+    CHARACTER          UPLO
+    INTEGER            INFO, LDA, N
+    REAL(dp)           A( LDA, * )
+    END SUBROUTINE
+
     REAL(dp) FUNCTION DLAMCH( CMACH )
     import :: dp
     CHARACTER          CMACH