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

Added mixed kernel with Woodbury 3x3 and update splitting. #25

Browse Source
This commit is contained in:
Francois Coppens 2021-07-26 12:19:29 +02:00
parent d2eb9f5a5d
commit decd977fff
1 changed files with 211 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

214
org/qmckl_sherman_morrison_woodbury.org
View File

@ -30,6 +30,7 @@ int main() {
#+end_src
* Sherman-Morrison Helper Functions
Helper functions that are used by the Sherman-Morrison-Woodbury
@ -1000,6 +1001,7 @@ rc = qmckl_sherman_morrison_splitting_c(context, splitting_Dim, splitting_N_upda
assert(rc == QMCKL_SUCCESS);
#+end_src
* Sherman-Morrison with Woodbury 2x2 and update splitting
This is like naïve Sherman-Morrising, but whenever a denominator is
@ -1095,7 +1097,7 @@ qmckl_exit_code qmckl_sherman_morrison_smw2s_c(const qmckl_context context,
uint64_t n_of_2blocks = N_updates / 2;
uint64_t remainder = N_updates % 2;
uint64_t length_2block = 2 * Dim;
uint64_t length_1block = 1 * Dim;
// uint64_t length_1block = 1 * Dim;
// Apply first 2*n_of_2blocks updates in n_of_2blocks blocks of 2 updates with
// Woodbury 2x2 kernel
@ -1204,6 +1206,211 @@ rc = qmckl_sherman_morrison_smw2s_c(context, smw2s_Dim, smw2s_N_updates,
assert(rc == QMCKL_SUCCESS);
#+end_src
* Sherman-Morrison with Woodbury 3x3 and 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.
** ~qmckl_sherman_morrison_smw3s~
:PROPERTIES:
:Name: qmckl_sherman_morrison_smw3s
:CRetType: qmckl_exit_code
: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.
#+NAME: qmckl_sherman_morrison_smw3s_args
| qmckl_context | context | in | Global state |
| uint64_t | Dim | in | Leading dimension of Slater_inv |
| uint64_t | N_updates | in | Number of rank-1 updates to be applied to Slater_inv |
| double | Updates[N_updates*Dim] | in | Array containing the updates |
| uint64_t | Updates_index[N_updates] | in | Array containing the rank-1 updates |
| double | Slater_inv[Dim*Dim] | inout | Array containing the inverse of a Slater-matrix |
*** Requirements
Add description of the input variables. (see for e.g. qmckl_distance.org)
*** C header
#+CALL: generate_c_header(table=qmckl_sherman_morrison_smw3s_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_smw3s_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_smw3s_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)
!logical, external :: qmckl_sherman_morrison_f
info = qmckl_sherman_morrison_smw3s(context, Dim, N_updates, Updates, Updates_index, Slater_inv)
end function qmckl_sherman_morrison_smw3s_f
#+end_src
*** Source C
#+begin_src c :tangle (eval c) :comments org
#include <stdbool.h>
#include "qmckl.h"
qmckl_exit_code qmckl_sherman_morrison_smw3s_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) {
// #ifdef DEBUG // Leave commented out since debugging information is not yet implemented in QMCkl.
// std::cerr << "Called qmckl_sherman_morrison_woodbury_3 with " << N_updates
// << " updates" << std::endl;
// #endif
qmckl_context local_context;
local_context = qmckl_context_create();
qmckl_exit_code rc;
uint64_t n_of_3blocks = N_updates / 3;
uint64_t remainder = N_updates % 3;
uint64_t length_3block = 3 * Dim;
// Apply first 3*n_of_3blocks updates in n_of_3blocks blocks of 3 updates with
// Woodbury 3x3 kernel
double later_updates[Dim * N_updates];
uint64_t later_index[N_updates];
uint64_t later = 0;
if (n_of_3blocks > 0) {
for (uint64_t i = 0; i < n_of_3blocks; i++) {
double *Updates_3block = &Updates[i * length_3block];
uint64_t *Updates_index_3block = &Updates_index[i * 3];
rc = qmckl_woodbury_3_c(local_context, Dim, Updates_3block, Updates_index_3block, Slater_inv);
if (rc != 0) { // Send the entire block to slagel_splitting
uint64_t l = 0;
rc = qmckl_slagel_splitting_c(Dim, 3, Updates_3block, Updates_index_3block,
Slater_inv, later_updates + (Dim * later), later_index + later, &l);
later = later + l;
}
}
}
if (remainder != 0) { // Apply last remaining block of 2 updates with Woodbury 2x2 kernel
double *Updates_remainder_block = &Updates[n_of_3blocks * length_3block];
uint64_t *Updates_index_remainder_block = &Updates_index[3 * n_of_3blocks];
uint64_t l = 0;
rc = qmckl_slagel_splitting_c(Dim, remainder, Updates_remainder_block, Updates_index_remainder_block,
Slater_inv, later_updates + (Dim * later), later_index + later, &l);
later = later + l;
}
if (later > 0) {
rc = qmckl_sherman_morrison_splitting_c(local_context, Dim, later, later_updates, later_index, Slater_inv);
}
}
#+end_src
*** Performance...
** C interface :noexport:
#+CALL: generate_c_interface(table=qmckl_sherman_morrison_smw3s_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_smw3s &
(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_smw3s_c
info = qmckl_sherman_morrison_smw3s_c &
(context, Dim, N_updates, Updates, Updates_index, Slater_inv)
end function qmckl_sherman_morrison_smw3s
#+end_src
#+CALL: generate_f_interface(table=qmckl_sherman_morrison_smw3s_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_smw3s &
(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_smw3s
end interface
#+end_src
*** Test :noexport:
[TODO: FMJC] Write tests for the Sherman-Morrison part.
#+begin_src c :tangle (eval c_test)
const uint64_t smw32s_Dim = 3;
const uint64_t smw32s_N_updates = 3;
const uint64_t smw32s_Updates_index[3] = {1, 1, 1};
const double smw32s_Updates[9] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
double smw32s_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_smw3s_c(context, smw32s_Dim, smw32s_N_updates,
smw32s_Updates, smw32s_Updates_index, smw32s_Slater_inv);
assert(rc == QMCKL_SUCCESS);
#+end_src
* Sherman-Morrison with Woodbury 3x3, 2x2 and update splitting
This is like naïve Sherman-Morrising, but whenever a denominator is
@ -1300,8 +1507,8 @@ qmckl_exit_code qmckl_sherman_morrison_smw32s_c(const qmckl_context context,
uint64_t n_of_3blocks = N_updates / 3;
uint64_t remainder = N_updates % 3;
uint64_t length_3block = 3 * Dim;
uint64_t length_2block = 2 * Dim;
uint64_t length_1block = 1 * Dim;
// uint64_t length_2block = 2 * Dim;
// uint64_t length_1block = 1 * Dim;
// Apply first 3*n_of_3blocks updates in n_of_3blocks blocks of 3 updates with
// Woodbury 3x3 kernel
@ -1422,6 +1629,7 @@ assert(rc == QMCKL_SUCCESS);
* End of files
#+begin_src c :comments link :tangle (eval c_test)