#+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 #include "qmckl.h" int main() { qmckl_context context; context = qmckl_context_create(); #+end_src * Sherman-Morrison Helper Functions [TODO: FMJC] Add doc ** ~qmckl_sherman_morrison_threshold~ :PROPERTIES: :Name: qmckl_sherman_morrison_threshold :CRetType: double :FRetType: double precision :END: [TODO: FMJC] Add doc #+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]; if (fabs(den) < threshold) { return false; } 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