TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-07-18 08:53:47 +02:00
Code Issues Projects Releases Wiki Activity

Splitting tests bug fix. #25

Browse Source
This commit is contained in:
vijay gopal chilkuri 2021-07-22 18:20:20 +02:00
parent e1b325ab18
commit 90b6333560
1 changed files with 104 additions and 25 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

99
org/qmckl_sherman_morrison_woodbury.org
View File

@ -126,8 +126,8 @@ return QMCKL_SUCCESS;
qmckl_exit_code qmckl_slagel_splitting_c(uint64_t Dim,
uint64_t N_updates,
double *Updates,
uint64_t *Updates_index,
const double *Updates,
const uint64_t *Updates_index,
double *Slater_inv,
double *later_updates,
uint64_t *later_index,
@ -213,7 +213,11 @@ qmckl_exit_code qmckl_slagel_splitting_c(uint64_t Dim,
:FRetType: qmckl_exit_code
:END:
This is the simplest of the available Sherman-Morrison-Woodbury kernels in QMCkl. It applies rank-1 updates one by one in the order that is given. It only checks if the denominator in the Sherman-Morrison formula is not too close to zero (and exit with an error if it does) during the application of an update.
This is the simplest of the available Sherman-Morrison-Woodbury
kernels in QMCkl. It applies rank-1 updates one by one in the order
that is given. It only checks if the denominator in the
Sherman-Morrison formula is not too close to zero (and exit with an
error if it does) during the application of an update.
#+NAME: qmckl_sherman_morrison_args
| qmckl_context | context | in | Global state |
@ -820,12 +824,16 @@ assert(rc == QMCKL_SUCCESS);
* Sherman-Morrison with update splitting
This is like naïve Sherman-Morrising, but whenever a denominator is found that is too close to zero the update
is split in half. Then one half is applied immediately and the other have is ket for later. When all the updates
have been processed, the list of split updates that have been kept for later are processed. If again applying an
update results in a denominator that is too close to zero, it is split in half again. One half is applied immediately
and one half is kept for later. The algorithm is done when no more updates have been kept for later. This recursion
will always end in a finite number of steps, unless the last original update causes a singular Slater-matrix.
This is like naïve Sherman-Morrising, but whenever a denominator is
found that is too close to zero the update is split in half. Then one
half is applied immediately and the other have is ket for later. When
all the updates have been processed, the list of split updates that
have been kept for later are processed. If again applying an update
results in a denominator that is too close to zero, it is split in
half again. One half is applied immediately and one half is kept for
later. The algorithm is done when no more updates have been kept for
later. This recursion will always end in a finite number of steps,
unless the last original update causes a singular Slater-matrix.
** ~qmckl_sherman_morrison_splitting~
:PROPERTIES:
@ -856,8 +864,34 @@ will always end in a finite number of steps, unless the last original update cau
#+CALL: generate_c_header(table=qmckl_sherman_morrison_splitting_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
#+RESULTS:
#+begin_src c :tangle (eval h_func) :comments org
qmckl_exit_code qmckl_sherman_morrison_splitting_c (
const qmckl_context context,
const uint64_t Dim,
const uint64_t N_updates,
const double* Updates,
const uint64_t* Updates_index,
double* Slater_inv );
#+end_src
*** Source Fortran
#+begin_src f90 :tangle (eval f)
integer function qmckl_sherman_morrison_splitting_f(context, Slater_inv, Dim, N_updates, &
Updates, Updates_index) result(info)
use qmckl
implicit none
integer(qmckl_context) , intent(in) :: context
integer*8 , intent(in), value :: Dim, N_updates
integer*8 , intent(in) :: Updates_index(N_updates)
real*8 , intent(in) :: Updates(N_updates*Dim)
real*8 , intent(inout) :: Slater_inv(Dim*Dim)
info = qmckl_sherman_morrison_splitting(context, Dim, N_updates, Updates, Updates_index, Slater_inv)
end function qmckl_sherman_morrison_splitting_f
#+end_src
*** Source C
#+begin_src c :tangle (eval c) :comments org
@ -898,8 +932,52 @@ qmckl_exit_code qmckl_sherman_morrison_splitting_c(const qmckl_context context,
#+CALL: generate_c_interface(table=qmckl_sherman_morrison_splitting_args,rettyp=get_value("FRetType"),fname=get_value("Name"))
#+RESULTS:
#+begin_src f90 :tangle (eval f) :comments org :exports none
integer(c_int32_t) function qmckl_sherman_morrison_splitting &
(context, Dim, N_updates, Updates, Updates_index, Slater_inv) &
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 :: Dim
integer (c_int64_t) , intent(in) , value :: N_updates
real (c_double ) , intent(in) :: Updates(N_updates*Dim)
integer (c_int64_t) , intent(in) :: Updates_index(N_updates)
real (c_double ) , intent(inout) :: Slater_inv(Dim*Dim)
integer(c_int32_t), external :: qmckl_sherman_morrison_splitting_c
info = qmckl_sherman_morrison_splitting_c &
(context, Dim, N_updates, Updates, Updates_index, Slater_inv)
end function qmckl_sherman_morrison_splitting
#+end_src
#+CALL: generate_f_interface(table=qmckl_sherman_morrison_splitting_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_sherman_morrison_splitting &
(context, Dim, N_updates, Updates, Updates_index, Slater_inv) &
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 :: Dim
integer (c_int64_t) , intent(in) , value :: N_updates
real (c_double ) , intent(in) :: Updates(N_updates*Dim)
integer (c_int64_t) , intent(in) :: Updates_index(N_updates)
real (c_double ) , intent(inout) :: Slater_inv(Dim*Dim)
end function qmckl_sherman_morrison_splitting
end interface
#+end_src
*** Test :noexport:
@ -908,13 +986,14 @@ qmckl_exit_code qmckl_sherman_morrison_splitting_c(const qmckl_context context,
#+begin_src c :tangle (eval c_test)
const uint64_t splitting_Dim = 3;
const uint64_t splitting_N_updates = 3;
const uint64_t splitting_Updates_index[3] = {1, 1, 1};
const double splitting_Updates[9] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
double splitting_Slater_inv[9] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
// [TODO : FMJC ] add realistic tests
rc = qmckl_sherman_morrison_splitting_c(context, splitting_Dim, splitting_Updates, splitting_Updates_index, splitting_Slater_inv);
rc = qmckl_sherman_morrison_splitting_c(context, splitting_Dim, splitting_N_updates, splitting_Updates, splitting_Updates_index, splitting_Slater_inv);
assert(rc == QMCKL_SUCCESS);
#+end_src