Skip to content
GitLab
Commit 1846b76b authored by Irina Lebedeva's avatar Irina Lebedeva
Browse files

Fixing use of SC and HC

parent 5cbcb62a
Hide whitespace changes
Inline Side-by-side
libOMM/src/omm.F90
View file @ 1846b76b
......@@ -112,7 +112,7 @@ subroutine omm(m,n,n_occ,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,scale_
type(matrix), save :: HWdd ! G^dag*H*G matrix (n x n)
type(matrix), allocatable, save :: SW(:) ! overlap matrix in WF basis (n x n)
! for each value of ip
type(matrix), save :: SC ! S*C matrix (n x m)
type(matrix), allocatable, save :: SC(:) ! S*C matrix (n x m) for each value of ip
type(matrix), save :: SG ! S*G matrix (n x m)
type(matrix), allocatable, save :: C_Chl(:) ! Cholesky-transformed WF coeffs.
! in orbital basis (n x m) for each value of ip
......@@ -215,7 +215,6 @@ subroutine omm(m,n,n_occ,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,scale_
if (HWdd%is_initialized) call m_deallocate(HWdd)
if (HWd%is_initialized) call m_deallocate(HWd)
if (SG%is_initialized) call m_deallocate(SG)
if (SC%is_initialized) call m_deallocate(SC)
if (allocated(P)) then
do i = 1, np
if (P(i)%is_initialized) call m_deallocate(P(i))
......@@ -228,6 +227,12 @@ subroutine omm(m,n,n_occ,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,scale_
end do
deallocate(C_Chl)
end if
if (allocated(SC)) then
do i = 1, np
if (SC(i)%is_initialized) call m_deallocate(SC(i))
end do
deallocate(SC)
end if
end if
if (.not. work2%is_initialized) call m_allocate(work2, n, m, label=m_storage,&
......@@ -282,6 +287,9 @@ subroutine omm(m,n,n_occ,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,scale_
! allocate matrices involved
if (.not. allocated(HW)) allocate(HW(np))
if (.not. allocated(SW)) allocate(SW(np))
if (.not. (use_Cholesky .or. no_S)) then
if (.not. allocated(SC)) allocate(SC(np))
end if
if (use_Cholesky) then
if (.not. allocated(C_Chl)) allocate(C_Chl(np))
end if
......@@ -298,8 +306,8 @@ subroutine omm(m,n,n_occ,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,scale_
blocksize1=bl_n, blocksize2=bl_n)
if (.not. SW(ip)%is_initialized) call m_allocate(SW(ip), n, n, label=m_storage,&
blocksize1=bl_n, blocksize2=bl_n)
if ((.not. use_Cholesky) .or. (.not. no_S)) then
if (.not. SC%is_initialized) call m_allocate(SC, n, m, label=m_storage,&
if (.not. (use_Cholesky .or. no_S)) then
if (.not. SC(ip)%is_initialized) call m_allocate(SC(ip), n, m, label=m_storage,&
blocksize1=bl_n, blocksize2=bl_m)
if (.not. SG%is_initialized) call m_allocate(SG, n, m, label=m_storage,&
blocksize1=bl_n, blocksize2=bl_m)
......@@ -346,12 +354,14 @@ subroutine omm(m,n,n_occ,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,scale_
call m_set(HC, 'a', 0.0_dp, 0.0_dp, m_operation)
call m_set(HG, 'a', 0.0_dp, 0.0_dp, m_operation)
call m_set(work1, 'a', 0.0_dp, 0.0_dp, m_operation)
if ((.not. use_Cholesky) .or. (.not. no_S)) then
if (.not. (use_Cholesky .or. no_S)) then
call m_set(SG, 'a', 0.0_dp, 0.0_dp, m_operation)
end if
call m_set(SW(ip), 'a', 0.0_dp, 0.0_dp, m_operation)
if ((.not. use_Cholesky) .or. (.not. no_S)) then
call m_set(SC, 'a', 0.0_dp, 0.0_dp, m_operation)
if(new_S .or. first_call(ip)) then
call m_set(SW(ip), 'a', 0.0_dp, 0.0_dp, m_operation)
if (.not. (use_Cholesky .or. no_S)) then
call m_set(SC(ip), 'a', 0.0_dp, 0.0_dp, m_operation)
end if
end if
if(sparse) then
......@@ -430,15 +440,16 @@ subroutine omm(m,n,n_occ,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,scale_
end if
! -calculate the overlap matrix in WF basis: SW=C^dag*S*C
! it is better to update it even if S has not changed
if (use_Cholesky) then
call mm_multiply(C_Chl(ip), 'n',&
C_Chl(ip), 'c', SW(ip), 1.0_dp, 0.0_dp, m_operation)
else if (no_S) then
call mm_multiply(C_min, 'n', C_min, 'c',&
SW(ip), 1.0_dp, 0.0_dp, m_operation)
else
call calc_AW(S, C_min, SW(ip), SC, m_operation)
if(new_S .or. first_call(ip)) then
if (use_Cholesky) then
call mm_multiply(C_Chl(ip), 'n',&
C_Chl(ip), 'c', SW(ip), 1.0_dp, 0.0_dp, m_operation)
else if (no_S) then
call mm_multiply(C_min, 'n', C_min, 'c',&
SW(ip), 1.0_dp, 0.0_dp, m_operation)
else
call calc_AW(S, C_min, SW(ip), SC(ip), m_operation)
end if
end if
! calculate the gradient: G=2*(2*H*C-S*C*HW-H*C*SW)
......@@ -448,7 +459,7 @@ subroutine omm(m,n,n_occ,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,scale_
else if (no_S) then
call calc_G(HW(ip), SW(ip), G, HC, C_min, m_operation, keep_sparsity=sparse)
else
call calc_G(HW(ip), SW(ip), G, HC, SC, m_operation, keep_sparsity=sparse)
call calc_G(HW(ip), SW(ip), G, HC, SC(ip), m_operation, keep_sparsity=sparse)
end if
! calculate the preconditioned gradient by premultiplying G by P
......@@ -463,7 +474,7 @@ subroutine omm(m,n,n_occ,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,scale_
! (the energy at C is given by the zeroth-order coeff. alpha(0))
if (use_precon) then
call mm_multiply(HC, 'n', PG, 'c', HWd, 1.0_dp, 0.0_dp, m_operation)
call mm_multiply(SC, 'n', PG, 'c', SWd, 1.0_dp, 0.0_dp, m_operation)
call mm_multiply(SC(ip), 'n', PG, 'c', SWd, 1.0_dp, 0.0_dp, m_operation)
else
call mm_multiply(HC, 'n', G, 'c', HWd, 1.0_dp, 0.0_dp, m_operation)
if (use_Cholesky) then
......@@ -471,7 +482,7 @@ subroutine omm(m,n,n_occ,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,scale_
else if (no_S) then
call mm_multiply(C_min, 'n', G, 'c', SWd, 1.0_dp, 0.0_dp, m_operation)
else
call mm_multiply(SC, 'n', G, 'c', SWd, 1.0_dp, 0.0_dp, m_operation)
call mm_multiply(SC(ip), 'n', G, 'c', SWd, 1.0_dp, 0.0_dp, m_operation)
end if
end if
......@@ -532,7 +543,7 @@ subroutine omm(m,n,n_occ,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,scale_
else if (no_S) then
call mm_multiply(C_min, 'n', D, 'c', SWd, 1.0_dp, 0.0_dp, m_operation)
else
call mm_multiply(SC, 'n', D, 'c', SWd, 1.0_dp, 0.0_dp, m_operation)
call mm_multiply(SC(ip), 'n', D, 'c', SWd, 1.0_dp, 0.0_dp, m_operation)
end if
call calc_AW(H, D, HWdd, HG, m_operation)
if (use_Cholesky .or. no_S) then
......@@ -578,6 +589,12 @@ subroutine omm(m,n,n_occ,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,scale_
call m_add(HWd, 'n', HW(ip), x_min, 1.0_dp, m_operation)
call m_add(HWdd, 'n', HW(ip), x_min**2, 1.0_dp, m_operation)
! recalculate HC ans SC
call m_add(HG, 'n', HC, x_min, 1.0_dp, m_operation)
if(.not. (use_Cholesky .or. no_S)) then
call m_add(SG, 'n', SC(ip), x_min, 1.0_dp, m_operation)
end if
e_diff = 2.0_dp * abs((e_min - e_min_old)/(e_min + e_min_old))
if ((mpi_rank == 0) .and. long_out) write(log_unit,&
'(a, 2(1x, i5), 2(1x, es15.7e3), 1x, a)') '|', i, j, e_min, e_diff, ' |'
......@@ -589,14 +606,12 @@ subroutine omm(m,n,n_occ,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,scale_
end if
! recalculate G at the minimum (or for the rescaled coeffs.)
call m_add(HG, 'n', HC, x_min, 1.0_dp, m_operation)
if (use_Cholesky) then
call calc_G(HW(ip), SW(ip), G, HC, C_Chl(ip), m_operation, keep_sparsity=sparse)
else if (no_S) then
call calc_G(HW(ip), SW(ip), G, HC, C_min, m_operation, keep_sparsity=sparse)
else
call m_add(SG, 'n', SC, x_min, 1.0_dp, m_operation)
call calc_G(HW(ip), SW(ip), G, HC, SC, m_operation, keep_sparsity=sparse)
call calc_G(HW(ip), SW(ip), G, HC, SC(ip), m_operation, keep_sparsity=sparse)
end if
if (use_precon) call mm_multiply(G, 'n', P(ip), 'n', PG, 1.0_dp, 0.0_dp, m_operation)
......@@ -634,7 +649,6 @@ subroutine omm(m,n,n_occ,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,scale_
if (G_p%is_initialized) call m_deallocate(G_p)
if (G%is_initialized) call m_deallocate(G)
if (SG%is_initialized) call m_deallocate(SG)
if (SC%is_initialized) call m_deallocate(SC)
if (SWdd%is_initialized) call m_deallocate(SWdd)
if (SWd%is_initialized) call m_deallocate(SWd)
if (HWdd%is_initialized) call m_deallocate(HWdd)
......
libOMM/src/omm_callback.F90
View file @ 1846b76b
......@@ -130,7 +130,7 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
type(matrix), save :: HWdd ! G^dag*H*G matrix (n x n)
type(matrix), allocatable, save :: SW(:) ! overlap matrix in WF basis (n x n)
! for each value of ip
type(matrix), save :: SC ! S*C matrix (n x m)
type(matrix), allocatable, save :: SC(:) ! S*C matrix (n x m) for each value of ip
type(matrix), save :: SG ! S*G matrix (n x m)
type(matrix), allocatable, save :: C_Chl(:) ! Cholesky-transformed WF coeffs.
! in orbital basis (n x m) for each value of ip
......@@ -233,7 +233,6 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
if (HWdd%is_initialized) call m_deallocate(HWdd)
if (HWd%is_initialized) call m_deallocate(HWd)
if (SG%is_initialized) call m_deallocate(SG)
if (SC%is_initialized) call m_deallocate(SC)
if (allocated(P)) then
do i = 1, np
if (P(i)%is_initialized) call m_deallocate(P(i))
......@@ -246,6 +245,12 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
end do
deallocate(C_Chl)
end if
if (allocated(SC)) then
do i = 1, np
if (SC(i)%is_initialized) call m_deallocate(SC(i))
end do
deallocate(SC)
end if
end if
if (.not. work2%is_initialized) call m_allocate(work2, n, m, label=m_storage,&
......@@ -300,6 +305,9 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
! allocate matrices involved
if (.not. allocated(HW)) allocate(HW(np))
if (.not. allocated(SW)) allocate(SW(np))
if (.not. (use_Cholesky .or. no_S)) then
if (.not. allocated(SC)) allocate(SC(np))
end if
if (use_Cholesky) then
if (.not. allocated(C_Chl)) allocate(C_Chl(np))
end if
......@@ -315,8 +323,8 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
blocksize1=bl_n, blocksize2=bl_n)
if (.not. SW(ip)%is_initialized) call m_allocate(SW(ip), n, n, label=m_storage,&
blocksize1=bl_n, blocksize2=bl_n)
if ((.not. use_Cholesky) .or. (.not. no_S)) then
if (.not. SC%is_initialized) call m_allocate(SC, n, m, label=m_storage,&
if (.not. (use_Cholesky .or. no_S)) then
if (.not. SC(ip)%is_initialized) call m_allocate(SC(ip), n, m, label=m_storage,&
blocksize1=bl_n, blocksize2=bl_m)
if (.not. SG%is_initialized) call m_allocate(SG, n, m, label=m_storage,&
blocksize1=bl_n, blocksize2=bl_m)
......@@ -363,12 +371,14 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
call m_set(HC, 'a', 0.0_dp, 0.0_dp, m_operation)
call m_set(HG, 'a', 0.0_dp, 0.0_dp, m_operation)
call m_set(work1, 'a', 0.0_dp, 0.0_dp, m_operation)
if ((.not. use_Cholesky) .or. (.not. no_S)) then
if (.not. (use_Cholesky .or. no_S)) then
call m_set(SG, 'a', 0.0_dp, 0.0_dp, m_operation)
end if
call m_set(SW(ip), 'a', 0.0_dp, 0.0_dp, m_operation)
if ((.not. use_Cholesky) .or. (.not. no_S)) then
call m_set(SC, 'a', 0.0_dp, 0.0_dp, m_operation)
if(new_S .or. first_call(ip)) then
call m_set(SW(ip), 'a', 0.0_dp, 0.0_dp, m_operation)
if (.not. (use_Cholesky .or. no_S)) then
call m_set(SC(ip), 'a', 0.0_dp, 0.0_dp, m_operation)
end if
end if
if(sparse) then
......@@ -448,7 +458,6 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
end if
! -calculate the overlap matrix in WF basis: SW=C^dag*S*C
! it is better to update it even if S has not changed
if (use_Cholesky) then
call mm_multiply(C_Chl(ip), 'n',&
C_Chl(ip), 'c', SW(ip), 1.0_dp, 0.0_dp, m_operation)
......@@ -456,7 +465,7 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
call mm_multiply(C_min, 'n', C_min, 'c',&
SW(ip), 1.0_dp, 0.0_dp, m_operation)
else
call calc_AW(S, C_min, SW(ip), SC, m_operation)
call calc_AW(S, C_min, SW(ip), SC(ip), m_operation)
end if
! calculate the gradient: G=2*(2*H*C-S*C*HW-H*C*SW)
......@@ -466,7 +475,7 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
else if (no_S) then
call calc_G(HW(ip), SW(ip), G, HC, C_min, m_operation, keep_sparsity=sparse)
else
call calc_G(HW(ip), SW(ip), G, HC, SC, m_operation, keep_sparsity=sparse)
call calc_G(HW(ip), SW(ip), G, HC, SC(ip), m_operation, keep_sparsity=sparse)
end if
! calculate the preconditioned gradient by premultiplying G by P
......@@ -481,7 +490,7 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
! (the energy at C is given by the zeroth-order coeff. alpha(0))
if (use_precon) then
call mm_multiply(HC, 'n', PG, 'c', HWd, 1.0_dp, 0.0_dp, m_operation)
call mm_multiply(SC, 'n', PG, 'c', SWd, 1.0_dp, 0.0_dp, m_operation)
call mm_multiply(SC(ip), 'n', PG, 'c', SWd, 1.0_dp, 0.0_dp, m_operation)
else
call mm_multiply(HC, 'n', G, 'c', HWd, 1.0_dp, 0.0_dp, m_operation)
if (use_Cholesky) then
......@@ -489,7 +498,7 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
else if (no_S) then
call mm_multiply(C_min, 'n', G, 'c', SWd, 1.0_dp, 0.0_dp, m_operation)
else
call mm_multiply(SC, 'n', G, 'c', SWd, 1.0_dp, 0.0_dp, m_operation)
call mm_multiply(SC(ip), 'n', G, 'c', SWd, 1.0_dp, 0.0_dp, m_operation)
end if
end if
......@@ -550,7 +559,7 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
else if (no_S) then
call mm_multiply(C_min, 'n', D, 'c', SWd, 1.0_dp, 0.0_dp, m_operation)
else
call mm_multiply(SC, 'n', D, 'c', SWd, 1.0_dp, 0.0_dp, m_operation)
call mm_multiply(SC(ip), 'n', D, 'c', SWd, 1.0_dp, 0.0_dp, m_operation)
end if
call calc_HW_callback(H, D, HWdd, HG, m_operation)
if (use_Cholesky .or. no_S) then
......@@ -596,6 +605,12 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
call m_add(HWd, 'n', HW(ip), x_min, 1.0_dp, m_operation)
call m_add(HWdd, 'n', HW(ip), x_min**2, 1.0_dp, m_operation)
! recalculate HC ans SC
call m_add(HG, 'n', HC, x_min, 1.0_dp, m_operation)
if(.not. (use_Cholesky .or. no_S)) then
call m_add(SG, 'n', SC(ip), x_min, 1.0_dp, m_operation)
end if
e_diff = 2.0_dp * abs((e_min - e_min_old)/(e_min + e_min_old))
if ((mpi_rank == 0) .and. long_out) write(log_unit,&
'(a, 2(1x, i5), 2(1x, es15.7e3), 1x, a)') '|', i, j, e_min, e_diff, ' |'
......@@ -607,14 +622,12 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
end if
! recalculate G at the minimum (or for the rescaled coeffs.)
call m_add(HG, 'n', HC, x_min, 1.0_dp, m_operation)
if (use_Cholesky) then
call calc_G(HW(ip), SW(ip), G, HC, C_Chl(ip), m_operation, keep_sparsity=sparse)
else if (no_S) then
call calc_G(HW(ip), SW(ip), G, HC, C_min, m_operation, keep_sparsity=sparse)
else
call m_add(SG, 'n', SC, x_min, 1.0_dp, m_operation)
call calc_G(HW(ip), SW(ip), G, HC, SC, m_operation, keep_sparsity=sparse)
call m_add(SG, 'n', SC(ip), x_min, 1.0_dp, m_operation)
end if
if (use_precon) call mm_multiply(G, 'n', P(ip), 'n', PG, 1.0_dp, 0.0_dp, m_operation)
......@@ -651,7 +664,6 @@ subroutine omm_callback(m,n,H,S,new_S,e_min,D_min,calc_ED,eta,C_min,init_C,T,sca
if (G_p%is_initialized) call m_deallocate(G_p)
if (G%is_initialized) call m_deallocate(G)
if (SG%is_initialized) call m_deallocate(SG)
if (SC%is_initialized) call m_deallocate(SC)
if (SWdd%is_initialized) call m_deallocate(SWdd)
if (SWd%is_initialized) call m_deallocate(SWd)
if (HWdd%is_initialized) call m_deallocate(HWdd)
......
Supports Markdown
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