#+TITLE: Sherman-Morrison-Woodbury #+SETUPFILE: ../tools/theme.setup #+INCLUDE: ../tools/lib.org Low- and high-level functions that use the Sherman-Morrison and Woodbury matrix inversion formulas to update the inverse of a non-singualr matrix * Headers #+begin_src elisp :noexport :results none :exports none (org-babel-lob-ingest "../tools/lib.org") #+end_src #+begin_src c :comments link :tangle (eval c_test) :noweb yes #include "qmckl.h" #include "assert.h" #ifdef HAVE_CONFIG_H #include "config.h" #endif #include int main() { qmckl_context context; context = qmckl_context_create(); qmckl_exit_code rc; #+end_src * Helper Functions Helper functions that are used by the Sherman-Morrison-Woodbury kernels. These functions can only be used internally by higher level functions. ** ~qmckl_slagel_splitting~ :PROPERTIES: :Name: qmckl_slagel_splitting :CRetType: double :FRetType: double precision :END: ~qmckl_slagel_splitting~ is used internally to apply a list of rank-1 updates while splitting an update if necessary. In case of a split it applies the first half of the update while putting the second half in waiting queue to be applied at the end. For a given update $u_j$ we define a threshold value $\epsilon_j$, which is the minimum value of $1+v_j^TS^{-1}u_j$ for a non-singular matrix $S$. If $1+v_j^TS^{-1}u_j \geq \epsilon_j$, the update is applied as usual. Otherwise, $u_j$ will be redefined as $\frac{u_j}{2}$, and the other half (to be applied at the end) will be defined using vectors $\frac{u_{j'}}{2}$ and $v_{j'}^T=v_{j'}^T$. #+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 | breakdown | in | Break-down parameter on which to fail or not | | 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 - ~Dim >= 2~ - ~N_updates >= 1~ - ~Updates~ is allocated with at least $1 \times 2 \times 8$ bytes - ~Updates_index~ is allocated with at least $1 \times 8$ bytes - ~breakdown~ is a small number such that $0 < breakdown << 1$ - ~Slater_inv~ is allocated with at least $Dim \times Dim \times 8$ bytes - ~later_updates~ is allocated with at least $1 \times 2 \times 8$ bytes - ~later_index~ is allocated with at least $1 \times 8$ bytes - ~later >= 0~ *** C header #+CALL: generate_c_header(table=qmckl_slagel_splitting_args,rettyp=get_value("CRetType"),fname=get_value("Name")) #+RESULTS: #+begin_src c :tangle (eval h_func) :comments org qmckl_exit_code qmckl_slagel_splitting_c ( 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 ); #+end_src *** Source Fortran *** Source C #+begin_src c :tangle (eval c) :comments org #include #include #include "qmckl.h" qmckl_exit_code qmckl_slagel_splitting_c(uint64_t Dim, 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) { // #ifdef DEBUG // Leave commented out since debugging information is not yet implemented in QMCkl. // std::cerr << "Called slagel_splitting with " << N_updates << " updates" << std::endl; // #endif double C[Dim]; double D[Dim]; uint64_t l = 0; // For each update while (l < N_updates) { // C = S^{-1} x U_l for (uint64_t i = 0; i < Dim; i++) { C[i] = 0; for (uint64_t j = 0; j < Dim; j++) { C[i] += Slater_inv[i * Dim + j] * Updates[l * Dim + j]; } } // Denominator double den = 1 + C[Updates_index[l] - 1]; if (fabs(den) < breakdown) { // U_l = U_l / 2 (do the split) for (uint64_t i = 0; i < Dim; i++) { later_updates[*later * Dim + i] = Updates[l * Dim + i] / 2.0; C[i] /= 2.0; } later_index[*later] = Updates_index[l]; (*later)++; den = 1 + C[Updates_index[l] - 1]; } double iden = 1 / den; // D = v^T x S^{-1} for (uint64_t j = 0; j < Dim; j++) { D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j]; } // S^{-1} = S^{-1} - C x D / den for (uint64_t i = 0; i < Dim; i++) { for (uint64_t j = 0; j < Dim; j++) { double update = C[i] * D[j] * iden; Slater_inv[i * Dim + j] -= update; } } l += 1; } return QMCKL_SUCCESS; } #+end_src *** Performance This function performce better for cycles with 1 rank-1 update and with a high fail-rate. ** C interface :noexport: #+CALL: generate_c_interface(table=qmckl_slagel_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_slagel_splitting & (Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, later_updates, later_index, later) & bind(C) result(info) use, intrinsic :: iso_c_binding implicit none 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(in) , value :: breakdown real (c_double ) , intent(inout) :: Slater_inv(Dim*Dim) real (c_double ) , intent(inout) :: later_updates(Dim * N_updates) integer (c_int64_t) , intent(inout) :: later_index(N_updates) integer (c_int64_t) , intent(inout) :: later integer(c_int32_t), external :: qmckl_slagel_splitting_c info = qmckl_slagel_splitting_c & (Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, later_updates, later_index, later) end function qmckl_slagel_splitting #+end_src *** Test :noexport: This function does not have an explicit test as it is only used internally by higher level Sherman-Morrison-Woodbury functions. * Naïve Sherman-Morrison ** ~qmckl_sherman_morrison~ :PROPERTIES: :Name: qmckl_sherman_morrison :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: This is the simplest of the available Sherman-Morrison-Woodbury kernels. 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 when an update is evaluated. It will exit with an error code of the denominator is too close to zero. The formula that is applied is \[ (S + uv^T)^{-1} = S^{-1} - \frac{S^{-1} uv^T S^{-1}}{1 + v^T S^{-1} u} \] where $S$ is the Slater-matrix, $u$ and $v^T$ are the column and row vectors containing the updates, $S^{-1}$ is the inverse of the Slater-matrix. #+NAME: qmckl_sherman_morrison_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 | breakdown | in | Break-down parameter on which to fail or not | | double | Slater_inv[Dim*Dim] | inout | Array containing the inverse of a Slater-matrix | *** Requirements - ~context~ is not ~QMCKL_NULL_CONTEXT~ - ~Dim >= 2~ - ~N_updates >= 1~ - ~Updates~ is allocated with at least $1 \times 2 \times 8$ bytes - ~Updates_index~ is allocated with at least $1 \times 8$ bytes - ~breakdown~ is a small number such that $0 < breakdown << 1$ - ~Slater_inv~ is allocated with at least $Dim \times Dim \times 8$ bytes *** C header #+CALL: generate_c_header(table=qmckl_sherman_morrison_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_c ( const qmckl_context context, const uint64_t Dim, const uint64_t N_updates, const double* Updates, const uint64_t* Updates_index, const double breakdown, double* Slater_inv ); #+end_src *** Source Fortran #+begin_src f90 :tangle (eval f) integer function qmckl_sherman_morrison_f(context, Dim, N_updates, & Updates, Updates_index, breakdown, Slater_inv) 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(in) :: breakdown real*8 , intent(inout) :: Slater_inv(Dim*Dim) !logical, external :: qmckl_sherman_morrison_f info = qmckl_sherman_morrison(context, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv) end function qmckl_sherman_morrison_f #+end_src *** Source C #+begin_src c :tangle (eval c) :comments org #include #include "qmckl.h" qmckl_exit_code qmckl_sherman_morrison_c(const qmckl_context context, const uint64_t Dim, const uint64_t N_updates, const double* Updates, const uint64_t* Updates_index, const double breakdown, double * Slater_inv) { // #ifdef DEBUG // std::cerr << "Called qmckl_sherman_morrison with " << N_updates << " updates" << std::endl; // #endif double C[Dim]; double D[Dim]; uint64_t l = 0; // For each update while (l < N_updates) { // C = A^{-1} x U_l for (uint64_t i = 0; i < Dim; i++) { C[i] = 0; for (uint64_t j = 0; j < Dim; j++) { C[i] += Slater_inv[i * Dim + j] * Updates[l * Dim + j]; } } // Denominator double den = 1 + C[Updates_index[l] - 1]; if (fabs(den) < breakdown) { return QMCKL_FAILURE; } double iden = 1 / den; // D = v^T x A^{-1} for (uint64_t j = 0; j < Dim; j++) { D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j]; } // A^{-1} = A^{-1} - C x D / den for (uint64_t i = 0; i < Dim; i++) { for (uint64_t j = 0; j < Dim; j++) { double update = C[i] * D[j] * iden; Slater_inv[i * Dim + j] -= update; } } l += 1; } return QMCKL_SUCCESS; } #+end_src *** Performance This function performs better when there is only 1 rank-1 update in the update cycle and the fail-rate of rank-1 updates is high. ** C interface :noexport: #+CALL: generate_c_interface(table=qmckl_sherman_morrison_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 & (context, Dim, N_updates, Updates, Updates_index, breakdown, 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(in) :: breakdown real (c_double ) , intent(inout) :: Slater_inv(Dim*Dim) integer(c_int32_t), external :: qmckl_sherman_morrison_c info = qmckl_sherman_morrison_c & (context, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv) end function qmckl_sherman_morrison #+end_src #+CALL: generate_f_interface(table=qmckl_sherman_morrison_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 & (context, Dim, N_updates, Updates, Updates_index, breakdown, 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(in) :: breakdown real (c_double ) , intent(inout) :: Slater_inv(Dim*Dim) end function qmckl_sherman_morrison 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 Dim = 2; const uint64_t N_updates = 2; const uint64_t Updates_index[2] = {0, 0}; const double Updates[4] = {0.0, 0.0, 0.0, 0.0}; const double breakdown = 1e-3; double Slater_inv[4] = {0.0, 0.0, 0.0, 0.0}; // [TODO : FMJC ] add realistic tests rc = qmckl_sherman_morrison_c(context, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv); assert(rc == QMCKL_SUCCESS); #+end_src * Woodbury 2x2 ** ~qmckl_woodbury_2~ :PROPERTIES: :Name: qmckl_woodbury_2 :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: The simplest version of the generalised Sherman-Morrison-Woodbury kernels. It is used to apply two rank-1 updates at once. The formula used in this algorithm is called the Woodbury Matrix Identity \[ (S + U V)^{-1} = S^{-1} - C B^{-1} D \] where $S$ is the Slater-matrix $U$ and $V$ are the matrices containing the updates and the canonical basis matrix $S^{-1}$ is the inverse of the Slater-matrix $C:= S^{-1}U$, a Dim $\times 2$ matrix $B := 1 + VC$, the $2 \times 2$ matrix that is going to be inverted $D := VS^{-1}$, a $2 \times Dim$ matrix #+NAME: qmckl_woodbury_2_args | qmckl_context | context | in | Global state | | uint64_t | Dim | in | Leading dimension of Slater_inv | | double | Updates[2*Dim] | in | Array containing the updates | | uint64_t | Updates_index[2] | in | Array containing the rank-1 updates | | double | breakdown | in | Break-down parameter on which to fail or not | | double | Slater_inv[Dim*Dim] | inout | Array containing the inverse of a Slater-matrix | *** Requirements - ~context~ is not ~qmckl_null_context~ - ~dim >= 2~ - ~updates~ is allocated with at least $2 \times 2 \times 8$ bytes - ~updates_index~ is allocated with $2 \times 8$ bytes - ~breakdown~ is a small number such that $0 < breakdown << 1$ - ~slater_inv~ is allocated with at least $dim \times dim \times 8$ bytes *** C header #+CALL: generate_c_header(table=qmckl_woodbury_2_args,rettyp=get_value("CRetType"),fname=get_value("Name")) #+RESULTS: #+begin_src c :tangle (eval h_func) :comments org qmckl_exit_code qmckl_woodbury_2_c ( const qmckl_context context, const uint64_t Dim, const double* Updates, const uint64_t* Updates_index, const double breakdown, double* Slater_inv ); #+end_src *** Source Fortran #+begin_src f90 :tangle (eval f) integer function qmckl_woodbury_2_f(context, Dim, & Updates, Updates_index, breakdown, Slater_inv) result(info) use qmckl implicit none integer(qmckl_context) , intent(in) :: context integer*8 , intent(in), value :: Dim integer*8 , intent(in) :: Updates_index(2) real*8 , intent(in) :: Updates(2*Dim) real*8 , intent(in) :: breakdown real*8 , intent(inout) :: Slater_inv(Dim*Dim) !logical, external :: qmckl_woodbury_2_f info = qmckl_woodbury_2(context, Dim, Updates, Updates_index, breakdown, Slater_inv) end function qmckl_woodbury_2_f #+end_src *** Source C #+begin_src c :tangle (eval c) :comments org #include #include "qmckl.h" qmckl_exit_code qmckl_woodbury_2_c(const qmckl_context context, const uint64_t Dim, const double* Updates, const uint64_t* Updates_index, const double breakdown, double * Slater_inv) { /* C := S^{-1} * U, dim x 2 B := 1 + V * C, 2 x 2 D := V * S^{-1}, 2 x dim */ // #ifdef DEBUG // Leave commented out since debugging information is not yet implemented in QMCkl. // std::cerr << "Called Woodbury 2x2 kernel" << std::endl; // #endif const uint64_t row1 = (Updates_index[0] - 1); const uint64_t row2 = (Updates_index[1] - 1); // Compute C = S_inv * U !! NON-STANDARD MATRIX MULTIPLICATION BECAUSE // OF LAYOUT OF 'Updates' !! double C[2 * Dim]; for (uint64_t i = 0; i < Dim; i++) { for (uint64_t j = 0; j < 2; j++) { C[i * 2 + j] = 0; for (uint64_t k = 0; k < Dim; k++) { C[i * 2 + j] += Slater_inv[i * Dim + k] * Updates[Dim * j + k]; } } } // Compute B = 1 + V * C const double B0 = C[row1 * 2] + 1; const double B1 = C[row1 * 2 + 1]; const double B2 = C[row2 * 2]; const double B3 = C[row2 * 2 + 1] + 1; // Check if determinant of inverted matrix is not zero double det = B0 * B3 - B1 * B2; if (fabs(det) < breakdown) { return QMCKL_FAILURE; } // Compute B^{-1} with explicit formula for 2x2 inversion double Binv[4], idet = 1.0 / det; Binv[0] = idet * B3; Binv[1] = -1.0 * idet * B1; Binv[2] = -1.0 * idet * B2; Binv[3] = idet * B0; // Compute tmp = B^{-1} x (V.S^{-1}) double tmp[2 * Dim]; for (uint64_t i = 0; i < 2; i++) { for (uint64_t j = 0; j < Dim; j++) { tmp[i * Dim + j] = Binv[i * 2] * Slater_inv[row1 * Dim + j]; tmp[i * Dim + j] += Binv[i * 2 + 1] * Slater_inv[row2 * Dim + j]; } } // Compute (S + U V)^{-1} = S^{-1} - C x tmp for (uint64_t i = 0; i < Dim; i++) { for (uint64_t j = 0; j < Dim; j++) { Slater_inv[i * Dim + j] -= C[i * 2] * tmp[j]; Slater_inv[i * Dim + j] -= C[i * 2 + 1] * tmp[Dim + j]; } } return QMCKL_SUCCESS; } #+end_src *** Performance This function is most efficient when used in cases where there are only 2 rank-1 updates. ** C interface :noexport: #+CALL: generate_c_interface(table=qmckl_woodbury_2_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_woodbury_2 & (context, Dim, Updates, Updates_index, breakdown, 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 real (c_double ) , intent(in) :: Updates(2*Dim) integer (c_int64_t) , intent(in) :: Updates_index(2) real (c_double ) , intent(in) :: breakdown real (c_double ) , intent(inout) :: Slater_inv(Dim*Dim) integer(c_int32_t), external :: qmckl_woodbury_2_c info = qmckl_woodbury_2_c & (context, Dim, Updates, Updates_index, breakdown, Slater_inv) end function qmckl_woodbury_2 #+end_src #+CALL: generate_f_interface(table=qmckl_woodbury_2_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_woodbury_2 & (context, Dim, Updates, Updates_index, breakdown, 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 real (c_double ) , intent(in) :: Updates(2*Dim) integer (c_int64_t) , intent(in) :: Updates_index(2) real (c_double ) , intent(in) :: breakdown real (c_double ) , intent(inout) :: Slater_inv(Dim*Dim) end function qmckl_woodbury_2 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 woodbury_Dim = 2; const uint64_t woodbury_Updates_index[2] = {1, 1}; const double woodbury_Updates[4] = {1.0, 1.0, 1.0, 1.0}; const double woodbury_breakdown = 1e-3; double woodbury_Slater_inv[4] = {1.0, 1.0, 1.0, 1.0}; // [TODO : FMJC ] add realistic tests rc = qmckl_woodbury_2_c(context, woodbury_Dim, woodbury_Updates, woodbury_Updates_index, woodbury_breakdown, woodbury_Slater_inv); assert(rc == QMCKL_SUCCESS); #+end_src * Woodbury 3x3 ** ~qmckl_woodbury_3~ :PROPERTIES: :Name: qmckl_woodbury_3 :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: The 3x3 version of the Woodbury 2x2 kernel. It is used to apply three rank-1 updates at once. The formula used in this kernel is the same as for Woodbury 2x2, except for the sizes of the following matrices: $C:= S^{-1}U$, a Dim $\times 3$ matrix $B := 1 + VC$, the $3 \times 3$ matrix that is going to be inverted $D := VS^{-1}$, a $3 \times Dim$ matrix #+NAME: qmckl_woodbury_3_args | qmckl_context | context | in | Global state | | uint64_t | Dim | in | Leading dimension of Slater_inv | | double | Updates[3*Dim] | in | Array containing the updates | | uint64_t | Updates_index[3] | in | Array containing the rank-1 updates | | double | breakdown | in | Break-down parameter on which to fail or not | | double | Slater_inv[Dim*Dim] | inout | Array containing the inverse of a Slater-matrix | *** Requirements - ~context~ is not ~qmckl_null_context~ - ~dim >= 2~ - ~updates~ is allocated with at least $3 \times 2 \times 8$ bytes - ~updates_index~ is allocated with $3 \times 8$ bytes - ~breakdown~ is a small number such that $0 < breakdown << 1$ - ~slater_inv~ is allocated with at least $dim \times dim \times 8$ bytes *** C header #+CALL: generate_c_header(table=qmckl_woodbury_3_args,rettyp=get_value("CRetType"),fname=get_value("Name")) #+RESULTS: #+begin_src c :tangle (eval h_func) :comments org qmckl_exit_code qmckl_woodbury_3_c ( const qmckl_context context, const uint64_t Dim, const double* Updates, const uint64_t* Updates_index, const double breakdown, double* Slater_inv ); #+end_src *** Source Fortran #+begin_src f90 :tangle (eval f) integer function qmckl_woodbury_3_f(context, Dim, & Updates, Updates_index, breakdown, Slater_inv) result(info) use qmckl implicit none integer(qmckl_context) , intent(in) :: context integer*8 , intent(in), value :: Dim integer*8 , intent(in) :: Updates_index(3) real*8 , intent(in) :: Updates(3*Dim) real*8 , intent(in) :: breakdown real*8 , intent(inout) :: Slater_inv(Dim*Dim) !logical, external :: qmckl_woodbury_3_f info = qmckl_woodbury_3(context, Dim, Updates, Updates_index, breakdown, Slater_inv) end function qmckl_woodbury_3_f #+end_src *** Source C #+begin_src c :tangle (eval c) :comments org #include #include "qmckl.h" qmckl_exit_code qmckl_woodbury_3_c(const qmckl_context context, const uint64_t Dim, const double* Updates, const uint64_t* Updates_index, const double breakdown, double * Slater_inv) { /* C := S^{-1} * U, dim x 3 B := 1 + V * C, 3 x 3 D := V * S^{-1}, 3 x dim ,*/ // #ifdef DEBUG // Leave commented out since debugging information is not yet implemented in QMCkl. // std::cerr << "Called Woodbury 3x3 kernel" << std::endl; // #endif const uint64_t row1 = (Updates_index[0] - 1); const uint64_t row2 = (Updates_index[1] - 1); const uint64_t row3 = (Updates_index[2] - 1); // Compute C = S_inv * U !! NON-STANDARD MATRIX MULTIPLICATION BECAUSE // OF LAYOUT OF 'Updates' !! double C[3 * Dim]; for (uint64_t i = 0; i < Dim; i++) { for (uint64_t j = 0; j < 3; j++) { C[i * 3 + j] = 0; for (uint64_t k = 0; k < Dim; k++) { C[i * 3 + j] += Slater_inv[i * Dim + k] * Updates[Dim * j + k]; } } } // Compute B = 1 + V.C const double B0 = C[row1 * 3] + 1; const double B1 = C[row1 * 3 + 1]; const double B2 = C[row1 * 3 + 2]; const double B3 = C[row2 * 3]; const double B4 = C[row2 * 3 + 1] + 1; const double B5 = C[row2 * 3 + 2]; const double B6 = C[row3 * 3]; const double B7 = C[row3 * 3 + 1]; const double B8 = C[row3 * 3 + 2] + 1; // Check if determinant of B is not too close to zero double det; det = B0 * (B4 * B8 - B5 * B7) - B1 * (B3 * B8 - B5 * B6) + B2 * (B3 * B7 - B4 * B6); if (fabs(det) < breakdown) { return QMCKL_FAILURE; } // Compute B^{-1} with explicit formula for 3x3 inversion double Binv[9], idet = 1.0 / det; Binv[0] = (B4 * B8 - B7 * B5) * idet; Binv[1] = -(B1 * B8 - B7 * B2) * idet; Binv[2] = (B1 * B5 - B4 * B2) * idet; Binv[3] = -(B3 * B8 - B6 * B5) * idet; Binv[4] = (B0 * B8 - B6 * B2) * idet; Binv[5] = -(B0 * B5 - B3 * B2) * idet; Binv[6] = (B3 * B7 - B6 * B4) * idet; Binv[7] = -(B0 * B7 - B6 * B1) * idet; Binv[8] = (B0 * B4 - B3 * B1) * idet; // Compute tmp = B^{-1} x (V.S^{-1}) double tmp[3 * Dim]; for (uint64_t i = 0; i < 3; i++) { for (uint64_t j = 0; j < Dim; j++) { tmp[i * Dim + j] = Binv[i * 3] * Slater_inv[row1 * Dim + j]; tmp[i * Dim + j] += Binv[i * 3 + 1] * Slater_inv[row2 * Dim + j]; tmp[i * Dim + j] += Binv[i * 3 + 2] * Slater_inv[row3 * Dim + j]; } } // Compute (S + U V)^{-1} = S^{-1} - C x tmp for (uint64_t i = 0; i < Dim; i++) { for (uint64_t j = 0; j < Dim; j++) { Slater_inv[i * Dim + j] -= C[i * 3] * tmp[j]; Slater_inv[i * Dim + j] -= C[i * 3 + 1] * tmp[Dim + j]; Slater_inv[i * Dim + j] -= C[i * 3 + 2] * tmp[2 * Dim + j]; } } return QMCKL_SUCCESS; } #+end_src *** Performance... This function is most efficient when used in cases where there are only 3 rank-1 updates. ** C interface :noexport: #+CALL: generate_c_interface(table=qmckl_woodbury_3_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_woodbury_3 & (context, Dim, Updates, Updates_index, breakdown, 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 real (c_double ) , intent(in) :: Updates(3*Dim) integer (c_int64_t) , intent(in) :: Updates_index(3) real (c_double ) , intent(in) , value :: breakdown real (c_double ) , intent(inout) :: Slater_inv(Dim*Dim) integer(c_int32_t), external :: qmckl_woodbury_3_c info = qmckl_woodbury_3_c & (context, Dim, Updates, Updates_index, breakdown, Slater_inv) end function qmckl_woodbury_3 #+end_src #+CALL: generate_f_interface(table=qmckl_woodbury_3_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_woodbury_3 & (context, Dim, Updates, Updates_index, breakdown, 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 real (c_double ) , intent(in) :: Updates(3*Dim) integer (c_int64_t) , intent(in) :: Updates_index(3) real (c_double ) , intent(in) , value :: breakdown real (c_double ) , intent(inout) :: Slater_inv(Dim*Dim) end function qmckl_woodbury_3 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 woodbury3_Dim = 3; const uint64_t woodbury3_Updates_index[3] = {1, 1, 1}; const double woodbury3_Updates[9] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0}; const double woodbury3_breakdown = 1e-3; double woodbury3_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_woodbury_3_c(context, woodbury3_Dim, woodbury3_Updates, woodbury3_Updates_index, woodbury3_breakdown, woodbury3_Slater_inv); assert(rc == QMCKL_SUCCESS); #+end_src * Sherman-Morrison with update splitting ** ~qmckl_sherman_morrison_splitting~ :PROPERTIES: :Name: qmckl_sherman_morrison_splitting :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: This is a variation on the 'Naive' Sherman-Morrison kernel. Whenever the denominator $1+v_j^T S^{-1} u_j$ in the Sherman-Morrison formula is deemed to be too close to zero, the update $u_j$ is split in half: $u_j \rightarrow \frac{1}{1} u_j$. One half is applied immediately --necessarily increasing the value of the denominator because of the split-- while the other halve is put in a queue that will be applied when all the remaining updates have been treated. The kernel is executed recursively until the queue is eiter empty and all updates are applied successfully, or the size of the queue equals the number of initial updates. In the last case the Slater-matrix that would have resulted from applying the updates is un-invertable and therefore the kernel exits with an exit code. #+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 | breakdown | in | Break-down parameter on which to fail or not | | double | Slater_inv[Dim*Dim] | inout | Array containing the inverse of a Slater-matrix | *** Requirements - ~context~ is not ~QMCKL_NULL_CONTEXT~ - ~Dim >= 2~ - ~N_updates >= 1~ - ~Updates~ is allocated with at least $1 \times 2 \times 8$ bytes - ~Updates_index~ is allocated with at least $1 \times 8$ bytes - ~breakdown~ is a small number such that $0 < breakdown << 1$ - ~Slater_inv~ is allocated with at least $Dim \times Dim \times 8$ bytes *** C header #+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, const double breakdown, double* Slater_inv ); #+end_src *** Source Fortran #+begin_src f90 :tangle (eval f) integer function qmckl_sherman_morrison_splitting_f(context, Dim, N_updates, & Updates, Updates_index, breakdown, Slater_inv) 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(in) :: breakdown real*8 , intent(inout) :: Slater_inv(Dim*Dim) info = qmckl_sherman_morrison_splitting(context, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv) end function qmckl_sherman_morrison_splitting_f #+end_src *** Source C #+begin_src c :tangle (eval c) :comments org #include #include "qmckl.h" 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, const double breakdown, double * Slater_inv) { // #ifdef DEBUG // Leave commented out since debugging information is not yet implemented in QMCkl. // std::cerr << "Called qmckl_sherman_morrison_splitting with " << N_updates << " updates" << std::endl; // #endif qmckl_exit_code rc; double later_updates[Dim * N_updates]; uint64_t later_index[N_updates]; uint64_t later = 0; rc = qmckl_slagel_splitting_c(Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, later_updates, later_index, &later); if (later > 0) { rc = qmckl_sherman_morrison_splitting_c(context, Dim, later, later_updates, later_index, breakdown, Slater_inv); } return QMCKL_SUCCESS; } #+end_src *** Performance... This kernel performs best when there are only 1 rank-1 update cycles and/or when the fail-rate is high. ** C interface :noexport: #+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, breakdown, 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(in) , value :: breakdown 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, breakdown, 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, breakdown, 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(in) , value :: breakdown real (c_double ) , intent(inout) :: Slater_inv(Dim*Dim) end function qmckl_sherman_morrison_splitting 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 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}; const double splitting_breakdown = 1e-3; 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_N_updates, splitting_Updates, splitting_Updates_index, splitting_breakdown, splitting_Slater_inv); assert(rc == QMCKL_SUCCESS); #+end_src * Woodbury 2x2 with Sherman-Morrison and update splitting ** ~qmckl_sherman_morrison_smw2s~ :PROPERTIES: :Name: qmckl_sherman_morrison_smw2s :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: The Woodbury 2x2 kernel with Sherman-Morrison and update splitting combines the low-level Woodbury 2x2 kernel and Sherman-Morrison with update splitting. For a given number of updates $N$ it splits the number of updates in blocks of two updates. The blocks of two updates are then applied one by one using Woodbury 2x2. If a block of updates fails, both updates in the block are applied with Sherman-Morrison instead, split if necessary and with their second half put in a queue. After all blocks are processed the remaining one update --in case there was an odd number of updates to begin with-- is also aplpied with Sherman-Morrison and split if necessary. The queue containing the collected second halves of all the processed updates is processed at the very end to avoid having many intermediate queues containing only a few updates that risks an increased probability of artificially created non-singular intermediate matrices, resulting from division up the total number of updates in blocks of three. #+NAME: qmckl_sherman_morrison_smw2s_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 | breakdown | in | Break-down parameter on which to fail or not | | double | Slater_inv[Dim*Dim] | inout | Array containing the inverse of a Slater-matrix | *** Requirements - ~context~ is not ~QMCKL_NULL_CONTEXT~ - ~Dim >= 2~ - ~N_updates >= 1~ - ~Updates~ is allocated with at least $1 \times 2 \times 8$ bytes - ~Updates_index~ is allocated with at least $1 \times 8$ bytes - ~breakdown~ is a small number such that $0 < breakdown << 1$ - ~Slater_inv~ is allocated with at least $Dim \times Dim \times 8$ bytes *** C header #+CALL: generate_c_header(table=qmckl_sherman_morrison_smw2s_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_smw2s_c ( const qmckl_context context, const uint64_t Dim, const uint64_t N_updates, const double* Updates, const uint64_t* Updates_index, const double breakdown, double* Slater_inv ); #+end_src *** Source Fortran #+begin_src f90 :tangle (eval f) integer function qmckl_sherman_morrison_smw2s_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_smw2s (context, Dim, N_updates, Updates, Updates_index, Slater_inv) end function qmckl_sherman_morrison_smw2s_f #+end_src *** Source C #+begin_src c :tangle (eval c) :comments org #include #include "qmckl.h" qmckl_exit_code qmckl_sherman_morrison_smw2s_c(const qmckl_context context, const uint64_t Dim, const uint64_t N_updates, const double* Updates, const uint64_t* Updates_index, const double breakdown, double * Slater_inv) { // #ifdef DEBUG // Leave commented out since debugging information is not yet implemented in QMCkl. // std::cerr << "Called qmckl_sherman_morrison_woodbury_2 with " << N_updates // << " updates" << std::endl; // #endif qmckl_exit_code rc; uint64_t n_of_2blocks = N_updates / 2; uint64_t remainder = N_updates % 2; uint64_t length_2block = 2 * Dim; // Apply first 2*n_of_2blocks updates in n_of_2blocks blocks of 2 updates with // Woodbury 2x2 kernel double later_updates[Dim * N_updates]; uint64_t later_index[N_updates]; uint64_t later = 0; if (n_of_2blocks > 0) { for (uint64_t i = 0; i < n_of_2blocks; i++) { const double *Updates_2block = &Updates[i * length_2block]; const uint64_t *Updates_index_2block = &Updates_index[i * 2]; rc = qmckl_woodbury_2_c(context, Dim, Updates_2block, Updates_index_2block, breakdown, Slater_inv); if (rc != 0) { // Send the entire block to slagel_splitting uint64_t l = 0; rc = qmckl_slagel_splitting_c(Dim, 2, Updates_2block, Updates_index_2block, breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l); later = later + l; } } } if (remainder == 1) { // Apply last remaining update with slagel_splitting const double *Updates_1block = &Updates[n_of_2blocks * length_2block]; const uint64_t *Updates_index_1block = &Updates_index[2 * n_of_2blocks]; uint64_t l = 0; rc = qmckl_slagel_splitting_c(Dim, 1, Updates_1block, Updates_index_1block, breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l); later = later + l; } if (later > 0) { rc = qmckl_sherman_morrison_splitting_c(context, Dim, later, later_updates, later_index, breakdown, Slater_inv); } return QMCKL_SUCCESS; } #+end_src *** Performance... This kernel performs best for the case of two rank-1 update and a low fail-rate. ** C interface :noexport: #+CALL: generate_c_interface(table=qmckl_sherman_morrison_smw2s_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_smw2s & (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_smw2s_c info = qmckl_sherman_morrison_smw2s_c & (context, Dim, N_updates, Updates, Updates_index, Slater_inv) end function qmckl_sherman_morrison_smw2s #+end_src #+CALL: generate_f_interface(table=qmckl_sherman_morrison_smw2s_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_smw2s & (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_smw2s 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 smw2s_Dim = 3; const uint64_t smw2s_N_updates = 3; const uint64_t smw2s_Updates_index[3] = {1, 1, 1}; const double smw2s_Updates[9] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0}; const double smw2s_breakdown = 1e-3; double smw2s_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_smw2s_c(context, smw2s_Dim, smw2s_N_updates, smw2s_Updates, smw2s_Updates_index, smw2s_breakdown, smw2s_Slater_inv); assert(rc == QMCKL_SUCCESS); #+end_src * Woodbury 3x3 with Sherman-Morrison and update splitting ** ~qmckl_sherman_morrison_smw3s~ :PROPERTIES: :Name: qmckl_sherman_morrison_smw3s :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: The Woodbury 3x3 kernel with Sherman-Morrison and update splitting combines the low-level Woodbury 3x3 kernel and Sherman-Morrison with update splitting. It works the same as Woodbury 2x2 with Sherman-Morrison and update splitting, except that the updates are divided in blocks of three rank-1 updates instead of blocks of two rank-1 updates. #+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 | breakdown | in | Break-down parameter on which to fail or not | | double | Slater_inv[Dim*Dim] | inout | Array containing the inverse of a Slater-matrix | *** Requirements - ~context~ is not ~QMCKL_NULL_CONTEXT~ - ~Dim >= 2~ - ~N_updates >= 1~ - ~Updates~ is allocated with at least $1 \times 2 \times 8$ bytes - ~Updates_index~ is allocated with at least $1 \times 8$ bytes - ~breakdown~ is a small number such that $0 < breakdown << 1$ - ~Slater_inv~ is allocated with at least $Dim \times Dim \times 8$ bytes *** 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, const double breakdown, 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, const double breakdown, 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_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++) { const double *Updates_3block = &Updates[i * length_3block]; const uint64_t *Updates_index_3block = &Updates_index[i * 3]; rc = qmckl_woodbury_3_c(context, Dim, Updates_3block, Updates_index_3block, breakdown, 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, breakdown, 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 const double *Updates_remainder_block = &Updates[n_of_3blocks * length_3block]; const 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, breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l); later = later + l; } if (later > 0) { rc = qmckl_sherman_morrison_splitting_c(context, Dim, later, later_updates, later_index, breakdown, Slater_inv); } return QMCKL_SUCCESS; } #+end_src *** Performance... This kernel performs best for the case of three rank-1 update and a low fail-rate. ** 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 smw3s_Dim = 3; const uint64_t smw3s_N_updates = 3; const uint64_t smw3s_Updates_index[3] = {1, 1, 1}; const double smw3s_Updates[9] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0}; const double smw3s_breakdown = 1e-3; double smw3s_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, smw3s_Dim, smw3s_N_updates, smw3s_Updates, smw3s_Updates_index, smw3s_breakdown, smw3s_Slater_inv); assert(rc == QMCKL_SUCCESS); #+end_src * Woodbury 3x3 and 2x2 with Sherman-Morrison and update splitting ** ~qmckl_sherman_morrison_smw32s~ :PROPERTIES: :Name: qmckl_sherman_morrison_smw32s :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: The Woodbury 3x3 and 2x2 kernel with Sherman-Morrison and update splitting combines the low-level Woodbury 3x3 kernel, the Woobury 2x2 kernel and Sherman-Morrison with update splitting. It works the almost the same as Woodbury 3x3 with Sherman-Morrison and update splitting, except that when there is a remainder of two rank-1 updates, it is first tried with Woodbury 2x2 instead of sending them all to Sherman-Morrison with update splitting. For example, in the case of 5 updates the updates are applied in 1 block of 3 updates end 1 block of 2 updates. #+NAME: qmckl_sherman_morrison_smw32s_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 | breakdown | in | Break-down parameter on which to fail or not | | double | Slater_inv[Dim*Dim] | inout | Array containing the inverse of a Slater-matrix | *** Requirements - ~context~ is not ~QMCKL_NULL_CONTEXT~ - ~Dim >= 2~ - ~N_updates >= 1~ - ~Updates~ is allocated with at least $1 \times 2 \times 8$ bytes - ~Updates_index~ is allocated with at least $1 \times 8$ bytes - ~breakdown~ is a small number such that $0 < breakdown << 1$ - ~Slater_inv~ is allocated with at least $Dim \times Dim \times 8$ bytes *** C header #+CALL: generate_c_header(table=qmckl_sherman_morrison_smw32s_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_smw32s_c ( const qmckl_context context, const uint64_t Dim, const uint64_t N_updates, const double* Updates, const uint64_t* Updates_index, const double breakdown, double* Slater_inv ); #+end_src *** Source Fortran #+begin_src f90 :tangle (eval f) integer function qmckl_sherman_morrison_smw32s_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_smw32s(context, Dim, N_updates, Updates, Updates_index, Slater_inv) end function qmckl_sherman_morrison_smw32s_f #+end_src *** Source C #+begin_src c :tangle (eval c) :comments org #include #include "qmckl.h" qmckl_exit_code qmckl_sherman_morrison_smw32s_c(const qmckl_context context, const uint64_t Dim, const uint64_t N_updates, const double* Updates, const uint64_t* Updates_index, const double breakdown, 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_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++) { const double *Updates_3block = &Updates[i * length_3block]; const uint64_t *Updates_index_3block = &Updates_index[i * 3]; rc = qmckl_woodbury_3_c(context, Dim, Updates_3block, Updates_index_3block, breakdown, 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, breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l); later = later + l; } } } if (remainder == 2) { // Apply last remaining block of 2 updates with Woodbury 2x2 kernel const double *Updates_2block = &Updates[n_of_3blocks * length_3block]; const uint64_t *Updates_index_2block = &Updates_index[3 * n_of_3blocks]; rc = qmckl_woodbury_2_c(context, Dim, Updates_2block, Updates_index_2block, breakdown, Slater_inv); if (rc != 0) { // Send the entire block to slagel_splitting uint64_t l = 0; rc = qmckl_slagel_splitting_c(Dim, 2, Updates_2block, Updates_index_2block, breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l); later = later + l; } } else if (remainder == 1) { // Apply last remaining update with slagel_splitting const double *Updates_1block = &Updates[n_of_3blocks * length_3block]; const uint64_t *Updates_index_1block = &Updates_index[3 * n_of_3blocks]; uint64_t l = 0; rc = qmckl_slagel_splitting_c(Dim, 1, Updates_1block, Updates_index_1block, breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l); later = later + l; } if (later > 0) { rc = qmckl_sherman_morrison_splitting_c(context, Dim, later, later_updates, later_index, breakdown, Slater_inv); } return QMCKL_SUCCESS; } #+end_src *** Performance... This kernel performs best when the number of rank-1 updates is larger than 3 and fail-rates are low. ** C interface :noexport: #+CALL: generate_c_interface(table=qmckl_sherman_morrison_smw32s_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_smw32s & (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_smw32s_c info = qmckl_sherman_morrison_smw32s_c & (context, Dim, N_updates, Updates, Updates_index, Slater_inv) end function qmckl_sherman_morrison_smw32s #+end_src #+CALL: generate_f_interface(table=qmckl_sherman_morrison_smw32s_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_smw32s & (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_smw32s 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}; const double smw32s_breakdown = 1e-3; 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_smw32s_c(context, smw32s_Dim, smw32s_N_updates, smw32s_Updates, smw32s_Updates_index, smw32s_breakdown, smw32s_Slater_inv); assert(rc == QMCKL_SUCCESS); #+end_src * End of files #+begin_src c :comments link :tangle (eval c_test) assert (qmckl_context_destroy(context) == QMCKL_SUCCESS); return 0; } #+end_src # -*- mode: org -*- # vim: syntax=c