Commit 9770c06a authored by Sebastian Ohlmann's avatar Sebastian Ohlmann

Merge branch 'fix_lalg_inverter' into 'develop'

lalg_inverter

See merge request !884
parents 9a5dd00c 530a16bf
Pipeline #147110432 passed with stage
in 0 seconds
......@@ -230,10 +230,10 @@ contains
select case(cv%method)
case(CURV_METHOD_UNIFORM)
Jac(1:sb%dim, 1:sb%dim) = sb%rlattice_primitive(1:sb%dim, 1:sb%dim)
jdet = lalg_determinant(sb%dim, Jac, invert = .false.)
jdet = lalg_determinant(sb%dim, Jac, preserve_mat = .false.)
case(CURV_METHOD_GYGI)
call curv_gygi_jacobian(sb, cv%gygi, x, dummy, Jac)
jdet = M_ONE/lalg_determinant(sb%dim, Jac, invert = .false.)
jdet = M_ONE/lalg_determinant(sb%dim, Jac, preserve_mat = .false.)
case(CURV_METHOD_BRIGGS)
call curv_briggs_jacobian_inv(sb, cv%briggs, chi, Jac)
jdet = M_ONE
......@@ -242,7 +242,7 @@ contains
end do
case(CURV_METHOD_MODINE)
call curv_modine_jacobian_inv(sb, cv%modine, chi, dummy, Jac)
jdet = M_ONE*lalg_determinant(sb%dim, Jac, invert = .false.)
jdet = M_ONE*lalg_determinant(sb%dim, Jac, preserve_mat = .false.)
end select
SAFE_DEALLOCATE_A(Jac)
......
......@@ -104,7 +104,7 @@ contains
call berry_phase_matrix(st, mesh, noccst, ik, ik, gvector, matrix)
if(noccst > 0) then
det = lalg_determinant(noccst, matrix(1:noccst, 1:noccst), invert = .false.) ** st%smear%el_per_state
det = lalg_determinant(noccst, matrix(1:noccst, 1:noccst), preserve_mat=.false.) ** st%smear%el_per_state
else
det = M_ONE
end if
......
......@@ -71,14 +71,12 @@ module lalg_adv_oct_m
module procedure deigensolve, zeigensolve
end interface lalg_eigensolve
!> Note that lalg_determinant and lalg_inverter are just wrappers
!! over the same routine.
interface lalg_determinant
module procedure ddeterminant, zdeterminant
end interface lalg_determinant
interface lalg_inverter
module procedure ddeterminant, zdeterminant
module procedure dinverter, zinverter
end interface lalg_inverter
interface lalg_sym_inverter
......@@ -202,7 +200,6 @@ contains
CMPLX, allocatable :: evectors(:, :), zevalues(:)
FLOAT, allocatable :: evalues(:)
CMPLX :: deter
integer :: ii
......@@ -239,7 +236,7 @@ contains
ex(1:nn, 1:nn) = evectors(1:nn, 1:nn)
deter = lalg_inverter(nn, evectors)
call lalg_inverter(nn, evectors)
do ii = 1, nn
evectors(1:nn, ii) = zevalues(1:nn)*evectors(1:nn, ii)
......@@ -281,7 +278,6 @@ contains
CMPLX, allocatable :: evectors(:, :), zevalues(:)
FLOAT, allocatable :: evalues(:)
CMPLX :: deter
integer :: ii
......@@ -322,7 +318,7 @@ contains
ex(1:nn, 1:nn) = evectors(1:nn, 1:nn)
deter = lalg_inverter(nn, evectors)
call lalg_inverter(nn, evectors)
do ii = 1, nn
evectors(1:nn, ii) = zevalues(1:nn)*evectors(1:nn, ii)
......
......@@ -599,23 +599,29 @@ end subroutine X(lowest_eigensolve)
! ---------------------------------------------------------
!> Invert a real symmetric or complex Hermitian square matrix a
R_TYPE function X(determinant)(n, a, invert) result(d)
R_TYPE function X(determinant)(n, a, preserve_mat) result(d)
integer, intent(in) :: n
R_TYPE, intent(inout) :: a(:,:) !< (n,n)
logical, intent(in), optional :: invert
R_TYPE, target, intent(inout) :: a(:,:) !< (n,n)
logical, intent(in) :: preserve_mat
integer :: info, i
integer, allocatable :: ipiv(:)
R_TYPE, allocatable :: work(:)
R_TYPE, pointer :: tmp_a(:,:)
! No PUSH_SUB, called too often
ASSERT(n > 0)
SAFE_ALLOCATE(work(1:n)) ! query?
SAFE_ALLOCATE(ipiv(1:n))
call lapack_getrf(n, n, a(1, 1), lead_dim(a), ipiv(1), info)
if(preserve_mat) then
SAFE_ALLOCATE(tmp_a(1:n, 1:n))
tmp_a(1:n, 1:n) = a(1:n, 1:n)
else
tmp_a => a
end if
call lapack_getrf(n, n, tmp_a(1, 1), lead_dim(tmp_a), ipiv(1), info)
if(info < 0) then
write(message(1), '(5a, i5)') 'In ', TOSTRING(X(determinant)), ', LAPACK ', TOSTRING(X(getrf)), ' returned info = ', info
call messages_fatal(1)
......@@ -624,35 +630,76 @@ R_TYPE function X(determinant)(n, a, invert) result(d)
d = M_ONE
do i = 1, n
if(ipiv(i) /= i) then
d = - d*a(i, i)
d = - d*tmp_a(i, i)
else
d = d*a(i, i)
d = d*tmp_a(i, i)
end if
end do
if(optional_default(invert, .true.)) then
call lapack_getri(n, a(1, 1), lead_dim(a), ipiv(1), work(1), n, info)
if(info /= 0) then
write(message(1), '(5a, i5)') 'In ', TOSTRING(X(determinant)), ', LAPACK ', TOSTRING(X(getri)), ' returned info = ', info
SAFE_DEALLOCATE_A(ipiv)
if(preserve_mat) then
SAFE_DEALLOCATE_P(tmp_a)
end if
end function X(determinant)
! ---------------------------------------------------------
!> Invert a real symmetric or complex Hermitian square matrix a
subroutine X(inverter)(n, a, det)
integer, intent(in) :: n
R_TYPE, intent(inout) :: a(:,:) !< (n,n)
R_TYPE, optional, intent(out) :: det
integer :: info, i
integer, allocatable :: ipiv(:)
R_TYPE, allocatable :: work(:)
! No PUSH_SUB, called too often
ASSERT(n > 0)
SAFE_ALLOCATE(work(1:n)) ! query?
SAFE_ALLOCATE(ipiv(1:n))
call lapack_getrf(n, n, a(1, 1), lead_dim(a), ipiv(1), info)
if(info < 0) then
write(message(1), '(5a, i5)') 'In ', TOSTRING(X(determinant)), ', LAPACK ', TOSTRING(X(getrf)), ' returned info = ', info
call messages_fatal(1)
end if
if(present(det)) then
det = M_ONE
do i = 1, n
if(ipiv(i) /= i) then
det = - det*a(i, i)
else
det = det*a(i, i)
end if
end do
end if
call lapack_getri(n, a(1, 1), lead_dim(a), ipiv(1), work(1), n, info)
if(info /= 0) then
write(message(1), '(5a, i5)') 'In ', TOSTRING(X(determinant)), ', LAPACK ', TOSTRING(X(getri)), ' returned info = ', info
! http://www.netlib.org/lapack/explore-3.1.1-html/zgetri.f.html
!* INFO (output) INTEGER
!* = 0: successful exit
!* < 0: if INFO = -i, the i-th argument had an illegal value
!* > 0: if INFO = i, U(i,i) is exactly zero; the matrix is
!* singular and its inverse could not be computed.
if(info < 0) then
write(message(2), '(a,i5,a)') 'Argument number ', -info, ' had an illegal value.'
else
write(message(2), '(a,i5,a)') 'Diagonal element ', info, ' of U is 0; matrix is singular.'
end if
call messages_fatal(2)
end if
if(info < 0) then
write(message(2), '(a,i5,a)') 'Argument number ', -info, ' had an illegal value.'
else
write(message(2), '(a,i5,a)') 'Diagonal element ', info, ' of U is 0; matrix is singular.'
end if
call messages_fatal(2)
end if
SAFE_DEALLOCATE_A(work)
SAFE_DEALLOCATE_A(ipiv)
end function X(determinant)
end subroutine X(inverter)
! ---------------------------------------------------------
......
......@@ -224,7 +224,7 @@
type(controlfunction_t), intent(inout) :: par
integer :: i, mm, nn, n, j, k
FLOAT :: t, det, w1, dt
FLOAT :: t, w1, dt
FLOAT, allocatable :: neigenvec(:, :), eigenvec(:, :), eigenval(:)
type(tdf_t), allocatable :: fnn(:)
......@@ -363,7 +363,7 @@
par%utransf = transpose(eigenvec)
par%utransfi = par%utransf
det = lalg_inverter(par%dim, par%utransfi)
call lalg_inverter(par%dim, par%utransfi)
SAFE_DEALLOCATE_A(eigenvec)
SAFE_DEALLOCATE_A(eigenval)
......@@ -448,7 +448,7 @@
end do
par%utransf = transpose(eigenvec)
par%utransfi = par%utransf
det = lalg_inverter(n-1, par%utransfi)
call lalg_inverter(n-1, par%utransfi)
SAFE_DEALLOCATE_A(eigenvec)
SAFE_DEALLOCATE_A(neigenvec)
......
......@@ -125,7 +125,6 @@
type(states_elec_t), intent(in) :: psi_in
type(states_elec_t), intent(inout) :: chi_out
CMPLX :: zdet
CMPLX, allocatable :: cI(:), dI(:), mat(:, :, :), mm(:, :, :, :), mk(:, :), lambda(:, :)
CMPLX, allocatable :: zpsi(:, :), zchi(:, :)
integer :: ik, ist, jst, ia, ib, n_pairs, nst, kpoints, jj, idim, ip
......@@ -154,7 +153,7 @@
dI(ia) = zstates_elec_mpdotp(namespace, gr%mesh, tg%est%st, psi_in, mat)
if(abs(dI(ia)) > CNST(1.0e-12)) then
do ik = 1, kpoints
zdet = lalg_inverter(nst, mm(1:nst, 1:nst, ik, ia))
call lalg_inverter(nst, mm(1:nst, 1:nst, ik, ia))
end do
end if
call zstates_elec_matrix_swap(mat, tg%est%pair(ia))
......
......@@ -144,7 +144,7 @@ subroutine X(broyden_extrapolation)(this, coeff, d1, d2, d3, vin, vnew, iter_use
FLOAT, parameter :: w0 = CNST(0.01), ww = M_FIVE
integer :: i, j, k, l
R_TYPE :: gamma, determinant
R_TYPE :: gamma
R_TYPE, allocatable :: beta(:, :), work(:)
PUSH_SUB(X(broyden_extrapolation))
......@@ -186,7 +186,7 @@ subroutine X(broyden_extrapolation)(this, coeff, d1, d2, d3, vin, vnew, iter_use
end do
! invert matrix beta
determinant = lalg_inverter(iter_used, beta)
call lalg_inverter(iter_used, beta)
do i = 1, iter_used
work(i) = M_ZERO
......@@ -468,7 +468,7 @@ subroutine X(pulay_extrapolation)(this, d2, d3, vin, vout, vnew, iter_used, f, d
R_TYPE, intent(out) :: vnew(:, :, :)
integer :: i, j, k, l
R_TYPE :: alpha, determinant
R_TYPE :: alpha
R_TYPE, allocatable :: a(:, :)
PUSH_SUB(X(pulay_extrapolation))
......@@ -496,7 +496,7 @@ subroutine X(pulay_extrapolation)(this, d2, d3, vin, vout, vnew, iter_used, f, d
return
end if
determinant = lalg_inverter(iter_used, a)
call lalg_inverter(iter_used, a)
! compute new vector
vnew = vin
......
......@@ -179,7 +179,7 @@ R_TYPE function X(states_elec_mpmatrixelement_g)(namespace, mesh, st1, st2, opst
end do
end do
det = lalg_determinant(i1+k1, cc, invert = .true.)
call lalg_inverter(i1+k1, cc, det)
cc = det * transpose(cc)
! And now, apply Lowdin`s formula.
......@@ -231,7 +231,7 @@ R_TYPE function X(states_elec_mpmatrixelement_g)(namespace, mesh, st1, st2, opst
end do
! Get the matrix of cofactors.
zz = lalg_determinant(i1, cc, invert = .true.)
call lalg_inverter(i1, cc, zz)
cc = zz * transpose(cc)
! And now, apply Lowdin`s formula.
......@@ -339,10 +339,11 @@ R_TYPE function X(states_elec_mpdotp_g)(namespace, mesh, st1, st2, mat) result(d
end do
end do
dotp = dotp * (lalg_determinant(i1+k1, bb, invert = .false.)) ** st1%d%kweights(ik)
dotp = dotp * (lalg_determinant(i1+k1, bb, preserve_mat = .true.)) ** st1%d%kweights(ik)
if(i1 > 0) then
dotp = dotp * (lalg_determinant(i1, bb(1:i1, 1:i1), invert = .false.)) ** st1%d%kweights(ik)
dotp = dotp * (lalg_determinant(i1, bb(1:i1, 1:i1), preserve_mat = .false.)) ** st1%d%kweights(ik)
end if
SAFE_DEALLOCATE_A(bb)
end do
case(SPIN_POLARIZED, SPINORS)
......@@ -368,7 +369,7 @@ R_TYPE function X(states_elec_mpdotp_g)(namespace, mesh, st1, st2, mat) result(d
end do
end do
dotp = dotp * lalg_determinant(i1, bb, invert = .false.) ** st1%d%kweights(ik)
dotp = dotp * lalg_determinant(i1, bb, preserve_mat = .false.) ** st1%d%kweights(ik)
SAFE_DEALLOCATE_A(bb)
end if
......
......@@ -511,7 +511,6 @@ contains
character(len=80) :: filename, tmp
integer :: iunit, ik, ist, ik2, ispin
FLOAT :: determinant
PUSH_SUB(kdotp_write_eff_mass)
......@@ -545,7 +544,7 @@ contains
tmp = int2str(ist)
write(iunit,'(a, a, a, f12.8, a, a)') 'State #', trim(tmp), ', Energy = ', &
units_from_atomic(units_out%energy, st%eigenval(ist, ik)), ' ', units_abbrev(units_out%energy)
determinant = lalg_inverter(gr%sb%periodic_dim, kdotp_vars%eff_mass_inv(:, :, ist, ik), invert = .true.)
call lalg_inverter(gr%sb%periodic_dim, kdotp_vars%eff_mass_inv(:, :, ist, ik))
call output_tensor(iunit, kdotp_vars%eff_mass_inv(:, :, ist, ik), gr%sb%periodic_dim, unit_one)
end do
......
......@@ -716,7 +716,6 @@ contains
integer, intent(in) :: nidx
integer, intent(in) :: idx(:)
R_TYPE :: det
R_TYPE, allocatable :: tmp1(:, :), tmp2(:, :), tmp3(:, :)
type(profile_t), save :: prof
......@@ -729,7 +728,7 @@ contains
call states_elec_blockt_mul(gr%mesh, st, constr_start, constr_start, &
constr, constr, tmp1, xpsi1 = all_constr, xpsi2 = all_constr)
det = lalg_inverter(nconstr, tmp1, invert = .true.)
call lalg_inverter(nconstr, tmp1)
call states_elec_blockt_mul(gr%mesh, st, constr_start, vs_start, &
constr, vs, tmp2, xpsi1 = all_constr, xpsi2 = idx(1:nidx))
call lalg_gemm(nconstr, nidx, nconstr, R_TOTYPE(M_ONE), tmp1, tmp2, R_TOTYPE(M_ZERO), tmp3)
......
......@@ -461,7 +461,7 @@ contains
! And now, perform the necessary transformation.
ip(1:3, 1:3) = kick%pol(1:3, 1:3)
dump = lalg_inverter(3, ip)
call lalg_inverter(3, ip)
do is = 1, nspin
do ie = 1, energy_steps
sigma(:, :, ie, is) = matmul( transpose(ip), matmul(sigmap(:, :, ie, is), ip) )
......
......@@ -607,7 +607,7 @@ contains
logical :: file_exists
integer :: i, j, nspin, time_steps, lmax, nfiles, k, add_lm, l, m, max_add_lm
integer, allocatable :: iunit(:)
FLOAT :: dt, lambda, det, dump, o0
FLOAT :: dt, lambda, dump, o0
type(unit_t) :: mp_unit
FLOAT, allocatable :: q(:), mu(:), qq(:, :), c(:)
character(len=20) :: filename
......@@ -733,7 +733,7 @@ contains
end do
end do
det = lalg_inverter(nfiles, qq, invert = .true.)
call lalg_inverter(nfiles, qq)
mu = matmul(qq, c)
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment