diff --git a/org/qmckl_sherman_morrison_woodbury.org b/org/qmckl_sherman_morrison_woodbury.org index 254a780..79d491b 100644 --- a/org/qmckl_sherman_morrison_woodbury.org +++ b/org/qmckl_sherman_morrison_woodbury.org @@ -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 +#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)