TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2025-01-03 10:06:09 +01:00
Code Issues Projects Releases Wiki Activity

Strange Fortran type error...

Browse Source
This commit is contained in:
Francois Coppens 2023-02-21 19:27:23 +01:00
parent 8ba882675e
commit 8216f682b3
1 changed files with 217 additions and 146 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

363
org/qmckl_sherman_morrison_woodbury.org
View File

@ -132,7 +132,7 @@ end subroutine convert
#+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(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, lds) :: Inverse
@ -146,7 +146,25 @@ subroutine copy_back(Inverse, s_inv, lds, 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
end subroutine copy_back end subroutine copy_back_inv
#+end_src
#+begin_src f90 :tangle (eval f) :comment org :exports none
subroutine copy_back_lu(Later_updates, later_upds, lds, nupdates)
implicit none
integer*8 , intent(in) :: lds, nupdates
real*8 , intent(in) , dimension(nupdates, lds) :: Later_updates
real*8 , intent(out) :: later_upds(nupdates * lds)
integer*8 :: i, j
! Copy updated inverse back to s_inv
do i = 1, nupdates
do j = 1, lds
later_upds((i - 1) * lds + j) = Later_updates(i, j)
end do
end do
end subroutine copy_back_lu
#+end_src #+end_src
#+begin_src f90 :tangle (eval f) #+begin_src f90 :tangle (eval f)
@ -539,20 +557,21 @@ return '\n'.join(result)
~qmckl_sm_naive~ is a generic function that contains decision making logic that calls the proper kernel based on the used library configuration (~--enable-doc~ and ~--enable-hpc~) and the passed array dimensions ~LDS~ and ~Dim~. ~qmckl_sm_naive~ is a generic function that contains decision making logic that calls the proper kernel based on the used library configuration (~--enable-doc~ and ~--enable-hpc~) and the passed array dimensions ~LDS~ and ~Dim~.
#+begin_src c :tangle (eval c) :comments org :noweb yes #+begin_src c :tangle (eval c) :comments org :noweb yes
qmckl_exit_code qmckl_sm_naive(const qmckl_context context, qmckl_exit_code qmckl_sm_naive(const qmckl_context context,
const uint64_t LDS, const uint64_t LDS,
const uint64_t Dim, const uint64_t Dim,
const uint64_t N_updates, const uint64_t N_updates,
const double* Updates, const double* Updates,
const uint64_t* Updates_index, const uint64_t* Updates_index,
const double breakdown, const double breakdown,
double* Slater_inv, double* Slater_inv,
double* determinant) { double* determinant) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return qmckl_failwith(context, return qmckl_failwith(
QMCKL_NULL_CONTEXT, context,
"qmckl_sm_naive", QMCKL_NULL_CONTEXT,
NULL); "qmckl_sm_naive",
NULL);
} }
#ifdef HAVE_HPC #ifdef HAVE_HPC
@ -562,26 +581,28 @@ qmckl_exit_code qmckl_sm_naive(const qmckl_context context,
} }
} }
else { // Updating smaller sub-matrix else { // Updating smaller sub-matrix
return qmckl_sm_naive_hpc(context, return qmckl_sm_naive_hpc(
LDS, context,
Dim, LDS,
N_updates, Dim,
Updates, N_updates,
Updates_index, Updates,
breakdown, Updates_index,
Slater_inv, breakdown,
determinant); Slater_inv,
determinant);
} }
#else #else
return qmckl_sm_naive_doc(context, return qmckl_sm_naive_doc(
LDS, context,
Dim, LDS,
N_updates, Dim,
Updates, N_updates,
Updates_index, Updates,
breakdown, Updates_index,
Slater_inv, breakdown,
determinant); Slater_inv,
determinant);
#endif #endif
return QMCKL_FAILURE; return QMCKL_FAILURE;
@ -621,6 +642,32 @@ interface
end interface end interface
#+end_src #+end_src
#+CALL: generate_f_interface(table=qmckl_sm_naive_args,rettyp=get_value("FRetType"),fname=get_value("Name"))
#+RESULTS:
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_sm_naive &
(context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: LDS
integer (c_int64_t) , intent(in) , value :: Dim
integer (c_int64_t) , intent(in) , value :: N_updates
real (c_double ) , intent(in) :: Updates(N_updates*LDS)
integer (c_int64_t) , intent(in) :: Updates_index(N_updates)
real (c_double ) , intent(in) , value :: breakdown
real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS)
real (c_double ) , intent(inout) :: determinant
end function qmckl_sm_naive
end interface
#+end_src
#+CALL: generate_f_interface(table=qmckl_sm_naive_args,rettyp=get_value("FRetType"),fname="qmckl_sm_naive_doc") #+CALL: generate_f_interface(table=qmckl_sm_naive_args,rettyp=get_value("FRetType"),fname="qmckl_sm_naive_doc")
#+RESULTS: #+RESULTS:
@ -647,32 +694,6 @@ interface
end interface end interface
#+end_src #+end_src
#+CALL: generate_f_interface(table=qmckl_sm_naive_args,rettyp=get_value("FRetType"),fname="qmckl_sm_naive")
#+RESULTS:
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_sm_naive &
(context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) &
bind(C)
use, intrinsic :: iso_c_binding
import
implicit none
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: LDS
integer (c_int64_t) , intent(in) , value :: Dim
integer (c_int64_t) , intent(in) , value :: N_updates
real (c_double ) , intent(in) :: Updates(N_updates*LDS)
integer (c_int64_t) , intent(in) :: Updates_index(N_updates)
real (c_double ) , intent(in) , value :: breakdown
real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS)
real (c_double ) , intent(inout) :: determinant
end function qmckl_sm_naive
end interface
#+end_src
*** Performance *** Performance
This function performs best when there is only 1 rank-1 update in the update cycle. It is This function performs best when there is only 1 rank-1 update in the update cycle. It is
not useful to use Sherman-Morrison with update splitting for these cycles since splitting not useful to use Sherman-Morrison with update splitting for these cycles since splitting
@ -803,27 +824,32 @@ integer function qmckl_sm_splitting_core_doc_f( &
updates_index, & updates_index, &
breakdown, & breakdown, &
s_inv, & s_inv, &
later_updates, & later_upds, &
later_index, & Later_index, &
later, & Later, &
determinant) result(info) determinant) result(info)
use qmckl use qmckl
implicit none implicit none
integer*8 , intent(in) :: context integer*8 , intent(in) :: context
integer*8 , intent(in) :: lds, dim integer*8 , intent(in) :: lds, dim
integer*8 , intent(in) :: nupdates integer*8 , intent(in) :: nupdates
integer*8 , intent(in) :: updates_index(nupdates) integer*8 , intent(in) :: updates_index(nupdates)
real*8 , intent(in) :: upds(nupdates * lds) real*8 , intent(in) :: upds(nupdates * lds)
real*8 , intent(in) :: breakdown real*8 , intent(in) :: breakdown
real*8 , intent(inout) :: s_inv(dim * lds) real*8 , intent(inout) :: s_inv(dim * lds)
real*8 , intent(inout) :: later_updates(nupdates * lds) real*8 , intent(inout) :: determinant
integer*8 , intent(inout) :: later_index(nupdates) integer*8 , intent(inout) :: Later
integer*8 , intent(inout) :: later integer*8 , intent(inout) :: Later_index(nupdates)
real*8 , intent(inout) :: determinant real*8 , intent(inout) :: later_upds(nupdates * lds)
real*8 , dimension(lds, nupdates) :: Updates real*8 , dimension(nupdates, lds) :: Updates
real*8 , dimension(nupdates, lds) :: Later_updates
real*8 , dimension(dim, lds) :: Inverse real*8 , dimension(dim, lds) :: Inverse
real*8 , dimension(dim) :: C
real*8 , dimension(lds) :: D
real*8 :: denominator, idenominator, update
integer*8 :: i, j, l, row
info = QMCKL_FAILURE info = QMCKL_FAILURE
@ -831,15 +857,65 @@ integer function qmckl_sm_splitting_core_doc_f( &
info = QMCKL_INVALID_CONTEXT info = QMCKL_INVALID_CONTEXT
return return
endif endif
! Convert 'upds' and 's_inv' into the more easily readable Fortran ! Convert 'upds' and 's_inv' into the more easily readable Fortran
! 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)
! YET TO BE IMPLEMENTED
! Copy updated inverse back to s_inv l = 1;
call copy_back(Inverse, s_inv, lds, dim) ! For each update do...
do while (l < nupdates + 1)
! Compute C = S^{-1}U(l)
do i = 1, dim
C(i) = 0
do j = 1, dim
C(i) = C(i) + Inverse(i, j) * Updates(j, l)
end do
end do
! Compute denominator = 1 + V(l)^TC
row = updates_index(l)
denominator = 1 + C(row)
! If denominator is too close to zero:
! - Split update in 2 before storing in Later_updates
! - Split previously computed vector C in 2
! - Recompute the denominator
if (abs(denominator) < breakdown) then
do i = 1, dim
Later_updates(i, l) = Updates(i, l) / 2
C(i) = C(i) / 2
end do
Later_index(Later) = updates_index(l)
Later = Later + 1
denominator = 1 + C(row)
end if
idenominator = 1 / denominator
! Update det(S)
determinant = determinant * denominator
! selecting column: v_l^T * S_inv
D = Inverse(row, :)
! A^{-1} = A^{-1} - C x D / denominator
do i = 1, dim
do j = 1, dim
update = C(i) * D(j) * idenominator
Inverse(i, j) = Inverse(i, j) - update
end do
end do
l = l + 1
end do
! Copy updated inverse and later updates
! back to s_inv and later_upds
call copy_back_inv(Inverse, s_inv, lds, dim)
call copy_back_lu(Later_Updates, later_upds, lds, nupdates)
info = QMCKL_SUCCESS info = QMCKL_SUCCESS
@ -855,57 +931,6 @@ for C users and in the module file 'qmckl_f.F90' for Fortran users.
#+CALL: generate_c_interface(table=qmckl_sm_splitting_core_args,rettyp=get_value("CRetType"),fname="qmckl_sm_splitting_core_doc") #+CALL: generate_c_interface(table=qmckl_sm_splitting_core_args,rettyp=get_value("CRetType"),fname="qmckl_sm_splitting_core_doc")
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_sm_splitting_core_doc &
(context, &
LDS, &
Dim, &
N_updates, &
Updates, &
Updates_index, &
breakdown, &
Slater_inv, &
later_updates, &
later_index, &
later, &
determinant) &
bind(C) result(info)
use, intrinsic :: iso_c_binding
implicit none
integer (c_int64_t) , intent(in) , value :: context
integer (c_int64_t) , intent(in) , value :: LDS
integer (c_int64_t) , intent(in) , value :: Dim
integer (c_int64_t) , intent(in) , value :: N_updates
real (c_double ) , intent(in) :: Updates(N_updates*LDS)
integer (c_int64_t) , intent(in) :: Updates_index(N_updates)
real (c_double ) , intent(in) , value :: breakdown
real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS)
real (c_double ) , intent(inout) :: later_updates(N_updates*LDS)
integer (c_int64_t) , intent(inout) :: later_index(N_updates)
integer (c_int64_t) , intent(inout) :: later
real (c_double ) , intent(inout) :: determinant
integer(c_int32_t), external :: qmckl_sm_splitting_core_doc_f
info = qmckl_sm_splitting_core_doc_f &
(context, &
LDS, &
Dim, &
N_updates, &
Updates, &
Updates_index, &
breakdown, &
Slater_inv, &
later_updates, &
later_index, &
later, &
determinant)
end function qmckl_sm_splitting_core_doc
#+end_src
*** C headers (exposed in qmckl.h) *** C headers (exposed in qmckl.h)
#+CALL: generate_c_header(table=qmckl_sm_splitting_core_args,rettyp=get_value("CRetType"),fname=get_value("Name")) #+CALL: generate_c_header(table=qmckl_sm_splitting_core_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
@ -926,10 +951,29 @@ qmckl_exit_code qmckl_sm_splitting_core (
double* determinant ); double* determinant );
#+end_src #+end_src
#+CALL: generate_c_header(table=qmckl_sm_splitting_core_args,rettyp=get_value("CRetType"),fname="qmckl_sm_splitting_core_hpc")
#+RESULTS:
#+begin_src c :tangle (eval h_private_func) :comments org
qmckl_exit_code qmckl_sm_splitting_core_hpc (
const qmckl_context context,
const uint64_t LDS,
const uint64_t Dim,
const uint64_t N_updates,
const double* Updates,
const uint64_t* Updates_index,
const double breakdown,
double* Slater_inv,
double* later_updates,
uint64_t* later_index,
uint64_t* later,
double* determinant );
#+end_src
#+CALL: generate_c_header(table=qmckl_sm_splitting_core_args,rettyp=get_value("CRetType"),fname="qmckl_sm_splitting_core_doc") #+CALL: generate_c_header(table=qmckl_sm_splitting_core_args,rettyp=get_value("CRetType"),fname="qmckl_sm_splitting_core_doc")
#+RESULTS: #+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org #+begin_src c :tangle (eval h_func) :no-expand comments org
qmckl_exit_code qmckl_sm_splitting_core_doc ( qmckl_exit_code qmckl_sm_splitting_core_doc (
const qmckl_context context, const qmckl_context context,
const uint64_t LDS, const uint64_t LDS,
@ -1386,7 +1430,8 @@ able to do numerically correct computations, it does not do it in the most effic
not be used in real workloads. not be used in real workloads.
#+begin_src f90 :tangle (eval f) #+begin_src f90 :tangle (eval f)
integer function qmckl_sm_splitting_doc_f(context, & integer recursive function qmckl_sm_splitting_doc_f( &
context, &
lds, dim, & lds, dim, &
nupdates, & nupdates, &
upds, & upds, &
@ -1397,28 +1442,54 @@ integer function qmckl_sm_splitting_doc_f(context, &
use qmckl use qmckl
implicit none implicit none
integer*8 , intent(in) :: context integer*8 , intent(in) :: context
integer*8 , intent(in) :: lds, dim integer*8 , intent(in) :: lds, dim
integer*8 , intent(in) :: nupdates integer*8 , intent(in) :: nupdates
integer*8 , intent(in) :: updates_index(nupdates) integer*8 , intent(in) :: updates_index(nupdates)
real*8 , intent(in) :: upds(nupdates * lds) real*8 , intent(in) :: upds(nupdates * lds)
real*8 , intent(in) :: breakdown real*8 , intent(in) :: breakdown
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 integer*8 :: Later
real*8 , dimension(dim, lds) :: Inverse integer*8 , dimension(nupdates) :: Later_index
real*8 , dimension(nupdates * lds) :: Later_updates
info = QMCKL_FAILURE info = QMCKL_FAILURE
! Convert 'upds' and 's_inv' into the more easily readable Fortran if (context == QMCKL_NULL_CONTEXT) then
! matrices 'Updates' and 'Inverse'. info = QMCKL_INVALID_CONTEXT
call convert(upds, s_inv, Updates, Inverse, nupdates, lds, dim) return
endif
! YET TO BE IMPLEMENTED Later = 0
Later_index = 0
Later_updates = 0
! Copy updated inverse back to s_inv info = qmckl_sm_splitting_core_doc_f( &
call copy_back(Inverse, s_inv, lds, dim) context, &
lds, dim, &
nupdates, &
upds, &
updates_index, &
breakdown, &
s_inv, &
Later_updates, &
Later_index, &
Later, &
determinant)
if (Later > 0) then
info = qmckl_sm_splitting_doc_f( &
context, &
lds, dim, &
Later, &
Later_updates, &
Later_index, &
breakdown, &
s_inv, &
determinant)
end if
info = QMCKL_SUCCESS info = QMCKL_SUCCESS