diff --git a/org/qmckl_sherman_morrison_woodbury.org b/org/qmckl_sherman_morrison_woodbury.org index 93dd8bd..b53c1ab 100644 --- a/org/qmckl_sherman_morrison_woodbury.org +++ b/org/qmckl_sherman_morrison_woodbury.org @@ -31,7 +31,7 @@ int main() { * Sherman-Morrison Helper Functions - + Helper functions that are used by the Sherman-Morrison-Woodbury kernels. These functions should only be used in the context of these kernels. @@ -97,14 +97,14 @@ return QMCKL_SUCCESS; This function is used to set the threshold value that is used in the kernels to determine if a matrix is invertable or not. In the Sherman-Morrison kernels this is determined by comparing the denominator in the Sherman-Morrison formula to the value set in threshold. If the value is smaller than the threshold value it means the matrix is not invertable. In the Woodbury kernels the threshold value is compared with the value of the determinant of the update matrix. #+NAME: qmckl_slagel_splitting_args - | 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 rank-1 updates | - | uint64_t | Updates_index[N_updates] | in | Array containing positions of the rank-1 updates | - | double | Slater_inv[Dim*Dim] | inout | Array containing the inverse Slater-matrix | - | double | later_updates[Dim * N_updates] | inout | Array containing the split updates for later | - | uint64_t | later_index[N_updates] | inout | Array containing the positions of the split updates for later | - | uint64_t | later | inout | Number of split updates for later | + | 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 rank-1 updates | + | uint64_t | Updates_index[N_updates] | in | Array containing positions of the rank-1 updates | + | double | Slater_inv[Dim*Dim] | inout | Array containing the inverse Slater-matrix | + | double | later_updates[Dim * N_updates] | inout | Array containing the split updates for later | + | uint64_t | later_index[N_updates] | inout | Array containing the positions of the split updates for later | + | uint64_t | later | inout | Number of split updates for later | *** Requirements @@ -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: @@ -841,12 +849,12 @@ will always end in a finite number of steps, unless the last original update cau error if it does) during the application of an update. #+NAME: qmckl_sherman_morrison_splitting_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 | + | 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 @@ -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