#+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 (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 #ifndef THRESHOLD #define THRESHOLD 1e-3 #endif int main() { qmckl_context context; context = qmckl_context_create(); #+end_src * 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. ** ~qmckl_sherman_morrison_threshold~ :PROPERTIES: :Name: qmckl_sherman_morrison_threshold :CRetType: double :FRetType: double precision :END: 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_sherman_morrison_threshold_args | double | thresh | out | Threshold | *** Requirements Add description of the input variables. (see for e.g. qmckl_distance.org) *** C header #+CALL: generate_c_header(table=qmckl_sherman_morrison_threshold_args,rettyp=get_value("CRetType"),fname=get_value("Name")) #+RESULTS: #+begin_src c :tangle (eval h_func) :comments org // Sherman-Morrison-Woodbury break-down threshold #ifndef THRESHOLD #define THRESHOLD 1e-3 #endif qmckl_exit_code qmckl_sherman_morrison_threshold ( double* const thresh ); #+end_src *** Source Fortran #+begin_src f90 :tangle (eval f) integer function qmckl_sherman_morrison_threshold_f(thresh) result(info) use qmckl implicit none real*8 , intent(inout) :: thresh !logical, external :: qmckl_sherman_morrison_f info = qmckl_sherman_morrison_threshold(thresh) end function qmckl_sherman_morrison_threshold_f #+end_src *** Source C #+begin_src c :tangle (eval c) :comments org #include #include #include "qmckl.h" // Sherman-Morrison-Woodbury break-down threshold static double qmckl_shreman_morrison_threshold(double* const threshold) { *threshold = THRESHOLD; // #ifdef DEBUG // std::cerr << "Break-down threshold set to: " << threshold << std::endl; // #endif } #+end_src *** Performance ** C interface :noexport: #+CALL: generate_c_interface(table=qmckl_sherman_morrison_threshold_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_threshold & (thresh) & bind(C) result(info) use, intrinsic :: iso_c_binding implicit none real (c_double ) , intent(out) :: thresh integer(c_int32_t), external :: qmckl_sherman_morrison_threshold_f info = qmckl_sherman_morrison_threshold_f & (thresh) end function qmckl_sherman_morrison_threshold #+end_src #+CALL: generate_f_interface(table=qmckl_sherman_morrison_threshold_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_threshold & (thresh) & bind(C) use, intrinsic :: iso_c_binding import implicit none real (c_double ) , intent(out) :: thresh end function qmckl_sherman_morrison_threshold end interface #+end_src *** Test :noexport: [TODO: FMJC] Write tests for the Sherman-Morrison part. * Naïve Sherman-Morrison ** ~qmckl_sherman_morrison~ :PROPERTIES: :Name: qmckl_sherman_morrison :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: The Sherman-Morrison formula \begin{align} S_k^{-1} &= (S_l + U_k)^-1 \\ &= S_l^{-1} - \frac{S_l^{-1}U_kS_l}{1+\underline{v}_k^tS_l^{-1}\underline{u}_k} \end{align} #+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 | 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_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, double* Slater_inv ); #+end_src *** Source Fortran #+begin_src f90 :tangle (eval f) integer function qmckl_sherman_morrison_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) :: 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(context, Dim, N_updates, Updates, Updates_index, 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, double * Slater_inv) { // #ifdef DEBUG // std::cerr << "Called qmckl_sherman_morrison with " << N_updates << " updates" << std::endl; // #endif double C[Dim]; double D[Dim]; double threshold = 0.0; qmckl_sherman_morrison_threshold(&threshold); unsigned int l = 0; // For each update while (l < N_updates) { // C = A^{-1} x U_l for (unsigned int i = 0; i < Dim; i++) { C[i] = 0; for (unsigned int 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]; double const thresh; qmckl_sherman_morrison_threshold(&thresh); if (fabs(den) < thresh) { return QMCKL_FAILURE; } double iden = 1 / den; // D = v^T x A^{-1} for (unsigned int j = 0; j < Dim; j++) { D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j]; } // A^{-1} = A^{-1} - C x D / den for (unsigned int i = 0; i < Dim; i++) { for (unsigned int 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 ** 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, 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_c info = qmckl_sherman_morrison_c & (context, Dim, N_updates, Updates, Updates_index, 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, 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 end interface #+end_src *** Test :noexport: [TODO: FMJC] Write tests for the Sherman-Morrison part. * 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