#+TITLE: Sherman-Morrison-Woodbury #+SETUPFILE: ../tools/theme.setup #+INCLUDE: ../tools/lib.org #+STARTUP: content Low- and high-level functions that use the Sherman-Morrison and Woodbury matrix inversion formulas to update the inverse of a non-singular 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 #+NAME:kernel_generator_range #+begin_src python :noweb yes :exports none range(2, 22) #+end_src * Naïve Sherman-Morrison :PROPERTIES: :Name: qmckl_sherman_morrison_naive :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 for any update $u_j$ (index $j$ is suppresed for clarity) 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. Even though the Slater-matrix $S$ with all updates applied at once is invertable, during the course of applying updates to the inverse Slater-matrix $S^{-1}$ one-by-one it can happen that one of the intermediate inverse matrices $S^{-1}$ becomes singular. Therefore a global threshold value $\epsilon$ is defined that is used to evaluate each individual update $u_j$ when it is applied. This value sets the lower bound for which the denominator $1+v_j^TS^{-1}u_j$ is considered to be too small and will most probably result in a singular matrix $S$, or at least in an inverse of $S$ of very poor numerical quality. Therefore, when $1+v_j^TS^{-1}u_j \geq \epsilon$, the update is applied as usual and the kernel exits with return code \texttt{QMCKL_SUCCESS}. If $1+v_j^TS^{-1}u_j \leq \epsilon$ the update is rejected and the kernel exits with return code \texttt{QMCKL_FAILURE}. If the determinant of the Slater-matrix is passed, it will be updated to the determinant resulting from applying the updates to the original matrix. #+NAME: qmckl_sherman_morrison_naive_args | Variable | Type | In/Out | Description | |-----------------+-------------------------+--------+------------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~LDS~ | ~uint64_t~ | in | Leading dimension of Slater_inv | | ~Dim~ | ~uint64_t~ | in | Dimension of Slater_inv | | ~N_updates~ | ~uint64_t~ | in | Number of rank-1 updates to be applied to Slater_inv | | ~Updates~ | ~double[N_updates*LDS]~ | in | Array containing the updates | | ~Updates_index~ | ~uint64_t[N_updates]~ | in | Array containing the rank-1 updates | | ~breakdown~ | ~double~ | in | Break-down parameter on which to fail or not | | ~Slater_inv~ | ~double[Dim*LDS]~ | inout | Array containing the inverse of a Slater-matrix | | ~determinant~ | ~double~ | inout | Determinant of the Slater-matrix | ** Pedagogical kernel source (in Fortran) The following source code written in Fortran is inteded to illustrate how the kernel works. Even though the kernel is able to do numerically correct computations, it does not do it in the most efficient way possible. It should therefore not be used in real workloads. #+begin_src f90 :tangle (eval f) integer function qmckl_sherman_morrison_naive_doc_f(context, & LDS, Dim, & N_updates, & Updates, & Updates_index, & breakdown, & Slater_inv, & determinant) result(info) use qmckl implicit none integer*8 , intent(in) :: context integer*8 , intent(in) :: LDS, Dim integer*8 , intent(in) :: N_updates integer*8 , intent(in) :: Updates_index(N_updates) real*8 , intent(in) :: Updates(N_updates*LDS) real*8 , intent(in) :: breakdown real*8 , intent(inout) :: Slater_inv(Dim*LDS) real*8 , intent(inout) :: determinant info = 0 if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif write(*,*) "Function 'qmckl_sherman_morrison_naive_doc_f' does nothing for now..." info = QMCKL_SUCCESS end function qmckl_sherman_morrison_naive_doc_f #+end_src *** C interface to the pedagogical kernel The following interface block in Fortran makes sure that the pedagogical kernel, written in Fortran, can be called from C using the ~ISO_C_BINDING~. #+CALL: generate_c_interface(table=qmckl_sherman_morrison_naive_args,rettyp=get_value("CRetType"),fname="qmckl_sherman_morrison_naive_doc") #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none integer(c_int32_t) function qmckl_sherman_morrison_naive_doc & (context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) & 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 :: LDS 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*LDS) 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*LDS) real (c_double ) , intent(inout) :: determinant integer(c_int32_t), external :: qmckl_sherman_morrison_naive_doc_f info = qmckl_sherman_morrison_naive_doc_f & (context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) end function qmckl_sherman_morrison_naive_doc #+end_src ** Requirements * ~context~ is not ~QMCKL_NULL_CONTEXT~ * ~LDS >= 2~ * ~Dim >= 2~ * ~N_updates >= 1~ * ~Updates~ is allocated with $N_updates \times Dim$ elements * ~Updates_index~ is allocated with $N_updates$ elements * ~breakdown~ is a small number such that $0 < breakdown << 1$ * ~Slater_inv~ is allocated with $Dim \times Dim$ elements * ~determinant > 0~ ** C headers (exposed in qmckl.h) #+CALL: generate_c_header(table=qmckl_sherman_morrison_naive_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_naive ( const qmckl_context context, const uint64_t LDS, const uint64_t Dim, const uint64_t N_updates, const double* Updates, const uint64_t* Updates_index, const double breakdown, double* Slater_inv, double* determinant ); #+end_src #+CALL: generate_c_header(table=qmckl_sherman_morrison_naive_args,rettyp=get_value("CRetType"),fname="qmckl_sherman_morrison_naive_hpc") #+RESULTS: #+begin_src c :tangle (eval h_func) :comments org qmckl_exit_code qmckl_sherman_morrison_naive_hpc ( const qmckl_context context, const uint64_t LDS, const uint64_t Dim, const uint64_t N_updates, const double* Updates, const uint64_t* Updates_index, const double breakdown, double* Slater_inv, double* determinant ); #+end_src #+CALL: generate_c_header(table=qmckl_sherman_morrison_naive_args,rettyp=get_value("CRetType"),fname="qmckl_sherman_morrison_naive_doc") #+RESULTS: #+begin_src c :tangle (eval h_func) :comments org qmckl_exit_code qmckl_sherman_morrison_naive_doc ( const qmckl_context context, const uint64_t LDS, const uint64_t Dim, const uint64_t N_updates, const double* Updates, const uint64_t* Updates_index, const double breakdown, double* Slater_inv, double* determinant ); #+end_src ** C sources Common includes and macros used by all the Sherman-Morrison-Woodbury kernels. #+begin_src c :tangle (eval c) :comments org #include #include #include "qmckl.h" #include "config.h" // Order important because // __GNUC__ also set in ICC, ICX and CLANG // __clang__ also set in ICX #if defined(__INTEL_COMPILER) #define IVDEP _Pragma("ivdep") #define ALIGNED _Pragma("vector aligned") #elif defined(__INTEL_LLVM_COMPILER) #define IVDEP _Pragma("ivdep") #define ALIGNED _Pragma("vector aligned") #elif defined(__clang__) #define IVDEP _Pragma("clang loop vectorize(enable)") #define ALIGNED #elif defined(__GNUC__) #define IVDEP _Pragma("GCC ivdep") #define ALIGNED #endif #+end_src ~qmckl_sherman_morrison_naive_hpc~ is a high performance variation of ~qmckl_sherman_morrison_naive~ written in C. It is used in cases when ~Dim~ is smaller than the leading dimension ~LDS~, irrespective of whetether ~LDS~ includes zero padding to benefit from SIMD instructions or not. Cases like this include situations where one wants to apply updates to a square submatrix of the full matrix. It takes advantage of memory aligned data and assumes no data dependencies inside the loops. The loops are fully vectorised whenever ~Dim~ is an integer multiple of ~SIMD_LEGTH~. #+begin_src c :tangle (eval c) :comments org qmckl_exit_code qmckl_sherman_morrison_naive_hpc( const qmckl_context context, const uint64_t LDS, const uint64_t Dim, const uint64_t N_updates, const double* __restrict Updates, const uint64_t* __restrict Updates_index, const double breakdown, double* __restrict Slater_inv, double* __restrict determinant) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith( context, QMCKL_NULL_CONTEXT, "qmckl_sherman_morrison_naive_hpc", NULL); } double __attribute__((aligned(8))) C[Dim]; double __attribute__((aligned(8))) D[LDS]; 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.0f; IVDEP ALIGNED for (uint64_t j = 0; j < Dim; j++) { C[i] += Slater_inv[i * LDS + j] * Updates[l * LDS + j]; } } // Denominator: v_l^T * C const int cui = Updates_index[l] - 1; double den = 1.0f + C[cui]; if (fabs(den) < breakdown) return QMCKL_FAILURE; double iden = 1.0f / den; // Update det(A) if (determinant) *determinant *= den; // selecting column: v_l^T * S_inv IVDEP ALIGNED for (uint64_t j = 0; j < Dim; j++) { D[j] = Slater_inv[cui * LDS + j]; } // A^{-1} = A^{-1} - C x D / den for (uint64_t i = 0; i < Dim; i++) { IVDEP ALIGNED for (uint64_t j = 0; j < Dim; j++) { const double update = C[i] * D[j] * iden; Slater_inv[i * LDS + j] -= update; } } l += 1; } return QMCKL_SUCCESS; } #+end_src #+NAME:naive_template_code #+begin_src c static inline qmckl_exit_code qmckl_sherman_morrison_naive_{Dim}( const qmckl_context context, const uint64_t N_updates, const double* __restrict Updates, const uint64_t* __restrict Updates_index, const double breakdown, double* __restrict Slater_inv, double* __restrict determinant) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith(context, QMCKL_NULL_CONTEXT, "qmckl_sherman_morrison_naive_{Dim}", NULL); } #define D{Dim}_P ((1+({Dim}-1)/SIMD_LENGTH)*SIMD_LENGTH) double __attribute__((aligned(8))) C[{Dim}]; double __attribute__((aligned(8))) D[D{Dim}_P]; 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; IVDEP ALIGNED for (uint64_t j = 0; j < D{Dim}_P; j++) { C[i] += Slater_inv[i * D{Dim}_P + j] * Updates[l * D{Dim}_P + j]; } } // Denominator const int cui = Updates_index[l] - 1; double den = 1.0f + C[cui]; if (fabs(den) < breakdown) { return QMCKL_FAILURE; } double iden = 1.0f / den; // Update det(A) if (determinant) *determinant *= den; // selecting column: D = v_l^T * S_inv IVDEP ALIGNED for (uint64_t j = 0; j < D{Dim}_P; j++) { D[j] = Slater_inv[cui * D{Dim}_P + j]; } // A^{-1} = A^{-1} - C x D / den for (uint64_t i = 0; i < {Dim}; i++) { IVDEP ALIGNED for (uint64_t j = 0; j < D{Dim}_P; j++) { double update = C[i] * D[j] * iden; Slater_inv[i * D{Dim}_P + j] -= update; } } l += 1; } return QMCKL_SUCCESS; } #+end_src #+NAME:naive_kernel_generator #+begin_src python :noweb yes :exports none text=""" <> """ result = [] for Dim in <>: Dim=str(Dim) result.append(text.replace("{Dim}",Dim) ) return '\n'.join(result) #+end_src #+NAME:naive_switch-case_generator #+begin_src python :noweb yes :exports none text=""" case {Dim}: return qmckl_sherman_morrison_naive_{Dim}(context, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant); """ result = [] for Dim in <>: Dim=str(Dim) result.append(text.replace("{Dim}",Dim) ) return '\n'.join(result) #+end_src #+begin_src c :tangle (eval c) :comments org :noweb yes <> qmckl_exit_code qmckl_sherman_morrison_naive(const qmckl_context context, const uint64_t LDS, const uint64_t Dim, const uint64_t N_updates, const double* Updates, const uint64_t* Updates_index, const double breakdown, double* Slater_inv, double* determinant) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith(context, QMCKL_NULL_CONTEXT, "qmckl_sherman_morrison_naive", NULL); } #ifdef HAVE_HPC if (LDS == (1+(Dim-1)/SIMD_LENGTH)*SIMD_LENGTH) { // Most cases switch (Dim) { <> } } else { // When SIMD_LENGTH > 1, called with LDS == Dim AND Dim != (1+(Dim-1)/SIMD_LENGTH)*SIMD_LENGTH) return qmckl_sherman_morrison_naive_hpc(context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant); } #else return qmckl_sherman_morrison_naive_doc(context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant); #endif return QMCKL_FAILURE; } #+end_src ** Fortran interfaces (exposed in qmckl_f.F90) :PROPERTIES: :Name: qmckl_sherman_morrison_naive :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+CALL: generate_f_interface(table=qmckl_sherman_morrison_naive_args,rettyp=get_value("FRetType"),fname="qmckl_sherman_morrison_naive_hpc") #+RESULTS: #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface integer(c_int32_t) function qmckl_sherman_morrison_naive_hpc & (context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) & 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 :: LDS 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*LDS) 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*LDS) real (c_double ) , intent(inout) :: determinant end function qmckl_sherman_morrison_naive_hpc end interface #+end_src #+CALL: generate_f_interface(table=qmckl_sherman_morrison_naive_args,rettyp=get_value("FRetType"),fname="qmckl_sherman_morrison_naive_doc") #+RESULTS: #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface integer(c_int32_t) function qmckl_sherman_morrison_naive_doc & (context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) & 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 :: LDS 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*LDS) 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*LDS) real (c_double ) , intent(inout) :: determinant end function qmckl_sherman_morrison_naive_doc end interface #+end_src #+CALL: generate_f_interface(table=qmckl_sherman_morrison_naive_args,rettyp=get_value("FRetType"),fname="qmckl_sherman_morrison_naive") #+RESULTS: #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface integer(c_int32_t) function qmckl_sherman_morrison_naive & (context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) & 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 :: LDS 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*LDS) 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*LDS) real (c_double ) , intent(inout) :: determinant end function qmckl_sherman_morrison_naive end interface #+end_src ** Tests The tests for the kernels are executed on datasets that are extracted from a run of QMC=Chem on Benzene (21 spin-up/21 spin down electrons) using 329 unique alpha determinants. The tests are run such that the kernels reject the computed inverse whenever the computed intermediate determinants or denominators are smaller than 1e-3. This is the default value in QMC=Chem. The tests will return QMCKL_SUCCESS whenever all the elements of the final matrix $R=S.S^-1 - 1$ are smaller than the given tolerance value of 1e-3, and will return QMCKL_FAILURE if the values are larger than this tolerance value. #+begin_src c :tangle (eval c_test) const uint64_t Dim = 21; const uint64_t LDS = (1+(Dim-1)/SIMD_LENGTH)*SIMD_LENGTH; const double breakdown = 1e-3; const double tolerance = 1e-3; double res[441]; #include "sm_test.h" assert(Updates1 != NULL); assert(Updates_index1 != NULL); assert(Slater_inv1 != NULL); // original determinant of Slater1 (before applying updates) double det = 3.407025646103221e-10; rc = qmckl_sherman_morrison_naive(context, LDS, Dim, N_updates1, Updates1, Updates_index1, breakdown, Slater_inv1, &det); // Check that the determinant is updated properly assert(fabs(det + 4.120398385068217e-10) < 1e-15); for (unsigned int i = 0; i < Dim; i++) { for (unsigned int j = 0; j < Dim; j++) { res[i * Dim + j] = 0; for (unsigned int k = 0; k < Dim; k++) { res[i * Dim + j] += Slater1[i * Dim + k] * Slater_inv1[k * LDS + j]; } } } rc = QMCKL_SUCCESS; for (unsigned int i = 0; i < Dim; i++) { for (unsigned int j = 0; j < Dim; j++) { if (i == j && fabs(res[i * Dim + j] - 1) > tolerance) { rc = QMCKL_FAILURE; } if (i != j && fabs(res[i * Dim + j]) > tolerance) { rc = QMCKL_FAILURE; } } } 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