From dcd9fa866d7d1f2dfa7c5abb6e987bf0f966cf11 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Fri, 20 Oct 2023 11:18:53 +0200 Subject: [PATCH 1/5] Add an interface to DSYGVX --- src/lapack.f90 | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/lapack.f90 b/src/lapack.f90 index 75fee7d..e2f64cc 100644 --- a/src/lapack.f90 +++ b/src/lapack.f90 @@ -77,6 +77,18 @@ SUBROUTINE DSYGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, & REAL(dp) A( LDA, * ), B( LDB, * ), W( * ), WORK( * ) END SUBROUTINE + SUBROUTINE DSYGVX( ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB, & + VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, & + LWORK, IWORK, IFAIL, INFO ) + import :: dp + CHARACTER JOBZ, RANGE, UPLO + INTEGER IL, INFO, ITYPE, IU, LDA, LDB, LDZ, LWORK, M, N + REAL(dp) ABSTOL, VL, VU + INTEGER IFAIL( * ), IWORK( * ) + REAL(dp) A( LDA, * ), B( LDB, * ), W( * ), WORK( * ), & + Z( LDZ, * ) + END SUBROUTINE + REAL(dp) FUNCTION DLAMCH( CMACH ) import :: dp CHARACTER CMACH From 44d32e9f44ddf915d5d4dfe854df6fe1c7825133 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Fri, 20 Oct 2023 11:29:17 +0200 Subject: [PATCH 2/5] Use dsygvx --- src/linalg.f90 | 26 ++++++++++++++++++++++---- 1 file changed, 22 insertions(+), 4 deletions(-) diff --git a/src/linalg.f90 b/src/linalg.f90 index a12ee3b..3920a71 100644 --- a/src/linalg.f90 +++ b/src/linalg.f90 @@ -1,6 +1,6 @@ module linalg use types, only: dp - use lapack, only: dsyevd, dsygvd, dgesv + use lapack, only: dsyevd, dsygvd, dgesv, dsygvx implicit none private public eigh, solve @@ -24,8 +24,10 @@ subroutine deigh_generalized(Am, Bm, lam, c) integer :: n ! lapack variables integer :: lwork, liwork, info - integer, allocatable :: iwork(:) - real(dp), allocatable :: Bmt(:,:), work(:) + integer, allocatable :: iwork(:), ifail(:) + real(dp), allocatable :: Bmt(:,:), work(:), Z(:,:) + integer :: il, iu, M + real(dp) :: abstol ! solve n = size(Am,1) @@ -35,8 +37,24 @@ subroutine deigh_generalized(Am, Bm, lam, c) lwork = 1 + 6*n + 2*n**2 liwork = 3 + 5*n allocate(Bmt(n,n), work(lwork), iwork(liwork)) + allocate(ifail(n)) c = Am; Bmt = Bm ! Bmt temporaries overwritten by dsygvd - call dsygvd(1,'V','L',n,c,n,Bmt,n,lam,work,lwork,iwork,liwork,info) + !call dsygvd(1,'V','L',n,c,n,Bmt,n,lam,work,lwork,iwork,liwork,info) + il = 1 + iu = 7 + M = iu-il+1 + allocate(z(n,M)) + abstol = 1e-8_dp + call dsygvx(1,'V','I','L',n,c,n,Bmt,n, & + 0._dp, 0._dp, 1, 7, abstol, M, lam, z, n, work, & + lwork, iwork, ifail, info) + c = 0 + c(:,:M) = z + !SUBROUTINE DSYGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, & + ! LWORK, IWORK, LIWORK, INFO ) + !SUBROUTINE DSYGVX( ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB, & + ! VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, & + ! LWORK, IWORK, IFAIL, INFO ) if (info /= 0) then print *, "dsygvd returned info =", info if (info < 0) then From aa4085120786b6f0e08e981d140cbf1696abd5c8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Fri, 20 Oct 2023 11:33:03 +0200 Subject: [PATCH 3/5] Simplify --- src/linalg.f90 | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/src/linalg.f90 b/src/linalg.f90 index 3920a71..f457d28 100644 --- a/src/linalg.f90 +++ b/src/linalg.f90 @@ -25,7 +25,7 @@ subroutine deigh_generalized(Am, Bm, lam, c) ! lapack variables integer :: lwork, liwork, info integer, allocatable :: iwork(:), ifail(:) - real(dp), allocatable :: Bmt(:,:), work(:), Z(:,:) + real(dp), allocatable :: Bmt(:,:), work(:), Z(:,:), Amt(:,:) integer :: il, iu, M real(dp) :: abstol @@ -37,19 +37,18 @@ subroutine deigh_generalized(Am, Bm, lam, c) lwork = 1 + 6*n + 2*n**2 liwork = 3 + 5*n allocate(Bmt(n,n), work(lwork), iwork(liwork)) + allocate(Amt(n,n)) allocate(ifail(n)) - c = Am; Bmt = Bm ! Bmt temporaries overwritten by dsygvd + Amt = Am; Bmt = Bm ! Bmt temporaries overwritten by dsygvd !call dsygvd(1,'V','L',n,c,n,Bmt,n,lam,work,lwork,iwork,liwork,info) il = 1 iu = 7 M = iu-il+1 allocate(z(n,M)) - abstol = 1e-8_dp - call dsygvx(1,'V','I','L',n,c,n,Bmt,n, & - 0._dp, 0._dp, 1, 7, abstol, M, lam, z, n, work, & + abstol = 1e-4_dp + call dsygvx(1,'V','I','L',n,Amt,n,Bmt,n, & + 0._dp, 0._dp, 1, 7, abstol, M, lam, c, n, work, & lwork, iwork, ifail, info) - c = 0 - c(:,:M) = z !SUBROUTINE DSYGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, & ! LWORK, IWORK, LIWORK, INFO ) !SUBROUTINE DSYGVX( ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB, & From 3284be8dbc187c301be7e10d71b06b521c219aa0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Fri, 20 Oct 2023 11:33:19 +0200 Subject: [PATCH 4/5] Hardwire accurate_eigensolver = .false. --- src/dirac.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dirac.f90 b/src/dirac.f90 index 0fc99c0..957aeb7 100644 --- a/src/dirac.f90 +++ b/src/dirac.f90 @@ -234,7 +234,7 @@ subroutine Ffunc(x, y, eng) real(dp) :: E_dirac_shift integer :: idx logical :: accurate_eigensolver - accurate_eigensolver = .true. + accurate_eigensolver = .false. iter = iter + 1 print *, "SCF iteration:", iter Vin = reshape(x, shape(Vin)) From a9ce4b23c6f2157a852bbdd99dbee6a2d86e6a7e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Fri, 20 Oct 2023 11:42:02 +0200 Subject: [PATCH 5/5] Do not copy matrices --- src/linalg.f90 | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/linalg.f90 b/src/linalg.f90 index f457d28..6ffa550 100644 --- a/src/linalg.f90 +++ b/src/linalg.f90 @@ -36,17 +36,15 @@ subroutine deigh_generalized(Am, Bm, lam, c) call assert_shape(c, [n, n], "eigh", "c") lwork = 1 + 6*n + 2*n**2 liwork = 3 + 5*n - allocate(Bmt(n,n), work(lwork), iwork(liwork)) - allocate(Amt(n,n)) + allocate(work(lwork), iwork(liwork)) allocate(ifail(n)) - Amt = Am; Bmt = Bm ! Bmt temporaries overwritten by dsygvd !call dsygvd(1,'V','L',n,c,n,Bmt,n,lam,work,lwork,iwork,liwork,info) il = 1 iu = 7 M = iu-il+1 allocate(z(n,M)) abstol = 1e-4_dp - call dsygvx(1,'V','I','L',n,Amt,n,Bmt,n, & + call dsygvx(1,'V','I','L',n,Am,n,Bm,n, & 0._dp, 0._dp, 1, 7, abstol, M, lam, c, n, work, & lwork, iwork, ifail, info) !SUBROUTINE DSYGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, &