4
1
mirror of https://github.com/pfloos/quack synced 2024-06-02 11:25:28 +02:00
quack/src/QuAcK/DIIS_extrapolation.f90

62 lines
1.1 KiB
Fortran
Raw Normal View History

2019-03-19 10:13:33 +01:00
subroutine DIIS_extrapolation(n_err,n_e,n_diis,error,e,error_in,e_inout)
! Perform DIIS extrapolation
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: n_err,n_e
double precision,intent(in) :: error_in(n_err),error(n_err,n_diis),e(n_e,n_diis)
! Local variables
double precision :: rcond
double precision,allocatable :: A(:,:),b(:),w(:)
! Output variables
integer,intent(inout) :: n_diis
double precision,intent(inout):: e_inout(n_e)
! Memory allocaiton
allocate(A(n_diis+1,n_diis+1),b(n_diis+1),w(n_diis+1))
! Update DIIS "history"
call prepend(n_err,n_diis,error,error_in)
call prepend(n_e,n_diis,e,e_inout)
! Build A matrix
A(1:n_diis,1:n_diis) = matmul(transpose(error),error)
A(1:n_diis,n_diis+1) = -1d0
A(n_diis+1,1:n_diis) = -1d0
A(n_diis+1,n_diis+1) = +0d0
! Build x matrix
b(1:n_diis) = +0d0
b(n_diis+1) = -1d0
! Solve linear system
call linear_solve(n_diis+1,A,b,w,rcond)
! Extrapolate
if(rcond > 1d-14) then
e_inout(:) = matmul(w(1:n_diis),transpose(e(:,1:n_diis)))
else
n_diis = 0
endif
end subroutine DIIS_extrapolation