TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2025-04-29 20:04:50 +02:00
Code Issues Projects Releases Wiki Activity

Made documentation vectors, matrices, conversion functions square for easier code comprehenisibility.

Browse Source
This commit is contained in:
Francois Coppens 2023-02-28 16:10:15 +01:00
parent 41745409a0
commit bb8377edc8
1 changed files with 40 additions and 42 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

80
org/qmckl_sherman_morrison_woodbury.org
View File

@ -110,21 +110,21 @@ subroutine convert(upds, s_inv, Updates, Inverse, nupdates, lds, dim)
integer*8 , intent(in) :: lds, dim, nupdates integer*8 , intent(in) :: lds, dim, nupdates
real*8 , intent(in) :: upds(nupdates * lds) real*8 , intent(in) :: upds(nupdates * lds)
real*8 , intent(in) :: s_inv(dim * lds) real*8 , intent(in) :: s_inv(dim * lds)
real*8 , intent(out) , dimension(lds, nupdates) :: Updates real*8 , intent(out) , dimension(dim, nupdates) :: Updates
real*8 , intent(out) , dimension(dim, lds) :: Inverse real*8 , intent(out) , dimension(dim, dim) :: Inverse
integer*8 :: i, j integer*8 :: i, j
! Construct Updates: lds x nupdates ! Construct Updates: lds x nupdates
do i = 1, nupdates do i = 1, nupdates
do j = 1, lds do j = 1, dim
Updates(j, i) = upds((i - 1) * lds + j) Updates(j, i) = upds((i - 1) * lds + j)
end do end do
end do end do
! Construct Inverse: dim x lds ! Construct Inverse: dim x lds
do i = 1, dim do i = 1, dim
do j = 1, lds do j = 1, dim
Inverse(i, j) = s_inv((i - 1) * lds + j) Inverse(i, j) = s_inv((i - 1) * lds + j)
end do end do
end do end do
@ -135,14 +135,14 @@ end subroutine convert
subroutine copy_back_inv(Inverse, s_inv, lds, dim) subroutine copy_back_inv(Inverse, s_inv, lds, dim)
implicit none implicit none
integer*8 , intent(in) :: lds, dim integer*8 , intent(in) :: lds, dim
real*8 , intent(in) , dimension(dim, lds) :: Inverse real*8 , intent(in) , dimension(dim, dim) :: Inverse
real*8 , intent(out) :: s_inv(dim * lds) real*8 , intent(out) :: s_inv(dim * lds)
integer*8 :: i, j integer*8 :: i, j
! Copy updated inverse back to s_inv ! Copy updated inverse back to s_inv
do i = 1, dim do i = 1, dim
do j = 1, lds do j = 1, dim
s_inv((i - 1) * lds + j) = Inverse(i, j) s_inv((i - 1) * lds + j) = Inverse(i, j)
end do end do
end do end do
@ -150,17 +150,17 @@ end subroutine copy_back_inv
#+end_src #+end_src
#+begin_src f90 :tangle (eval f) :comment org :exports none #+begin_src f90 :tangle (eval f) :comment org :exports none
subroutine copy_back_lu(Later_updates, later_upds, lds, nupdates) subroutine copy_back_lu(Later_updates, later_upds, lds, dim, nupdates)
implicit none implicit none
integer*8 , intent(in) :: lds, nupdates integer*8 , intent(in) :: lds, dim, nupdates
real*8 , intent(in) , dimension(lds, nupdates) :: Later_updates real*8 , intent(in) , dimension(dim, nupdates) :: Later_updates
real*8 , intent(out) :: later_upds(nupdates * lds) real*8 , intent(out) :: later_upds(nupdates * lds)
integer*8 :: i, j integer*8 :: i, j
! Copy updated inverse back to s_inv ! Copy updated inverse back to s_inv
do i = 1, nupdates do i = 1, nupdates
do j = 1, lds do j = 1, dim
later_upds((i - 1) * lds + j) = Later_updates(j, i) later_upds((i - 1) * lds + j) = Later_updates(j, i)
end do end do
end do end do
@ -188,10 +188,10 @@ integer function qmckl_sm_naive_doc_f(context, &
real*8 , intent(inout) :: s_inv(dim * lds) real*8 , intent(inout) :: s_inv(dim * lds)
real*8 , intent(inout) :: determinant real*8 , intent(inout) :: determinant
real*8 , dimension(lds, nupdates) :: Updates real*8 , dimension(dim, nupdates) :: Updates
real*8 , dimension(dim, lds) :: Inverse real*8 , dimension(dim, dim) :: Inverse
real*8 , dimension(dim) :: C real*8 , dimension(dim) :: C
real*8 , dimension(lds) :: D real*8 , dimension(dim) :: D
real*8 :: denominator, idenominator, update real*8 :: denominator, idenominator, update
integer*8 :: i, j, l, row integer*8 :: i, j, l, row
@ -843,16 +843,14 @@ integer function qmckl_sm_splitting_core_doc_f( &
integer*8 , intent(inout) :: Later_index(nupdates) integer*8 , intent(inout) :: Later_index(nupdates)
real*8 , intent(inout) :: later_upds(lds * nupdates) real*8 , intent(inout) :: later_upds(lds * nupdates)
real*8 , dimension(lds, nupdates) :: Updates real*8 , dimension(dim, nupdates) :: Updates
real*8 , dimension(lds, nupdates) :: Later_updates real*8 , dimension(dim, nupdates) :: Later_updates
real*8 , dimension(dim, lds) :: Inverse real*8 , dimension(dim, dim) :: Inverse
real*8 , dimension(dim) :: C real*8 , dimension(dim) :: C
real*8 , dimension(lds) :: D real*8 , dimension(dim) :: D
real*8 :: denominator, idenominator, update real*8 :: denominator, idenominator, update
integer*8 :: i, j, l, row integer*8 :: i, j, l, row
write(*,*) "Entering 'qmckl_sm_splittinig_core_doc_f'"
info = QMCKL_FAILURE info = QMCKL_FAILURE
if (context == QMCKL_NULL_CONTEXT) then if (context == QMCKL_NULL_CONTEXT) then
@ -864,7 +862,6 @@ integer function qmckl_sm_splitting_core_doc_f( &
! matrices 'Updates' and 'Inverse'. ! matrices 'Updates' and 'Inverse'.
call convert(upds, s_inv, Updates, Inverse, nupdates, lds, dim) call convert(upds, s_inv, Updates, Inverse, nupdates, lds, dim)
l = 1; l = 1;
! For each update do... ! For each update do...
do while (l < nupdates + 1) do while (l < nupdates + 1)
@ -917,12 +914,10 @@ integer function qmckl_sm_splitting_core_doc_f( &
! Copy updated inverse and later updates ! Copy updated inverse and later updates
! back to s_inv and later_upds ! back to s_inv and later_upds
call copy_back_inv(Inverse, s_inv, lds, dim) call copy_back_inv(Inverse, s_inv, lds, dim)
call copy_back_lu(Later_Updates, later_upds, lds, nupdates) call copy_back_lu(Later_Updates, later_upds, lds, dim, nupdates)
info = QMCKL_SUCCESS info = QMCKL_SUCCESS
write(*,*) "Leaving 'qmckl_sm_splittinig_core_doc_f'"
end function qmckl_sm_splitting_core_doc_f end function qmckl_sm_splitting_core_doc_f
#+end_src #+end_src
@ -1492,8 +1487,8 @@ integer function qmckl_woodbury_2x2_doc_f(&
real*8 , dimension(lds, 2) :: Updates real*8 , dimension(lds, 2) :: Updates
real*8 , dimension(dim, lds) :: Inverse real*8 , dimension(dim, lds) :: Inverse
real*8 , dimension(dim) :: C real*8 , dimension(dim, 2) :: C
real*8 , dimension(lds) :: D real*8 , dimension(2, dim) :: D
real*8 :: denominator, idenominator, update real*8 :: denominator, idenominator, update
integer*8 :: i, j, l, row integer*8 :: i, j, l, row
@ -1508,6 +1503,13 @@ integer function qmckl_woodbury_2x2_doc_f(&
! matrices 'Updates' and 'Inverse'. ! matrices 'Updates' and 'Inverse'.
call convert(upds, s_inv, Updates, Inverse, int(2,8), lds, dim) call convert(upds, s_inv, Updates, Inverse, int(2,8), lds, dim)
! Compute C(dim,2) = Inverse(dim,dim) x Updates(dim,2)
do i = 1, dim
do j = 1, dim
C(i,1) = C(i,1) + Inverse(1,j) * Updates(j,1)
C(i,2) = C(i,1) + Inverse(1,j) * Updates(j,2)
end do
end do
! Copy updated inverse and later updates ! Copy updated inverse and later updates
@ -1828,7 +1830,16 @@ qmckl_exit_code qmckl_woodbury_2x2(const qmckl_context context,
determinant); determinant);
} }
#else #else
return qmckl_woodbury_2x2_doc( // return qmckl_woodbury_2x2_doc(
// context,
// LDS,
// Dim,
// Updates,
// Updates_index,
// breakdown,
// Slater_inv,
// determinant);
return qmckl_woodbury_2x2_hpc(
context, context,
LDS, LDS,
Dim, Dim,
@ -1988,9 +1999,7 @@ integer recursive function qmckl_sm_splitting_doc_f( &
integer*8 :: Later integer*8 :: Later
integer*8 , dimension(nupdates) :: Later_index integer*8 , dimension(nupdates) :: Later_index
real*8 , dimension(lds * nupdates) :: Later_updates real*8 , dimension(dim * nupdates) :: Later_updates
write(*,*) "Entering 'qmckl_sm_splitting_doc_f'"
info = QMCKL_FAILURE info = QMCKL_FAILURE
@ -2030,8 +2039,6 @@ integer recursive function qmckl_sm_splitting_doc_f( &
info = QMCKL_SUCCESS info = QMCKL_SUCCESS
write(*,*) "Leaving 'qmckl_sm_splitting_doc_f'"
end function qmckl_sm_splitting_doc_f end function qmckl_sm_splitting_doc_f
#+end_src #+end_src
@ -2057,21 +2064,16 @@ integer(c_int32_t) function qmckl_sm_splitting_doc &
integer (c_int64_t) , intent(in) , value :: LDS integer (c_int64_t) , intent(in) , value :: LDS
integer (c_int64_t) , intent(in) , value :: Dim integer (c_int64_t) , intent(in) , value :: Dim
integer (c_int64_t) , intent(in) , value :: N_updates integer (c_int64_t) , intent(in) , value :: N_updates
real (c_double ) , intent(in) :: Updates(LDS*N_updates) real (c_double ) , intent(in) :: Updates(N_updates*LDS)
integer (c_int64_t) , intent(in) :: Updates_index(N_updates) integer (c_int64_t) , intent(in) :: Updates_index(N_updates)
real (c_double ) , intent(in) , value :: breakdown real (c_double ) , intent(in) , value :: breakdown
real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS) real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS)
real (c_double ) , intent(inout) :: determinant real (c_double ) , intent(inout) :: determinant
integer(c_int32_t), external :: qmckl_sm_splitting_doc_f integer(c_int32_t), external :: qmckl_sm_splitting_doc_f
write(*,*) "Entering 'qmckl_sm_splitting_doc'"
info = qmckl_sm_splitting_doc_f & info = qmckl_sm_splitting_doc_f &
(context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) (context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant)
write(*,*) "Leaving 'qmckl_sm_splitting_doc'"
end function qmckl_sm_splitting_doc end function qmckl_sm_splitting_doc
#+end_src #+end_src
@ -2215,8 +2217,6 @@ qmckl_exit_code qmckl_sm_splitting(
double* Slater_inv, double* Slater_inv,
double* determinant) { double* determinant) {
printf("Entering 'qmckl_sm_splitting'\n");
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith( return qmckl_failwith(
context, context,
@ -2248,8 +2248,6 @@ qmckl_exit_code qmckl_sm_splitting(
determinant); determinant);
#endif #endif
printf("Leaving 'qmckl_sm_splitting'\n");
return QMCKL_SUCCESS; return QMCKL_SUCCESS;
} }
#+end_src #+end_src