From 94e9b1396369b27f382c708429046cd48fc41af3 Mon Sep 17 00:00:00 2001 From: Pablo Oliveira Date: Wed, 13 Oct 2021 10:45:54 +0200 Subject: [PATCH] Update determinant in Woodbury 2x2 and fix tests --- org/qmckl_sherman_morrison_woodbury.org | 30 ++++++++++++++++++------- 1 file changed, 22 insertions(+), 8 deletions(-) diff --git a/org/qmckl_sherman_morrison_woodbury.org b/org/qmckl_sherman_morrison_woodbury.org index eaf4477..87a27d8 100644 --- a/org/qmckl_sherman_morrison_woodbury.org +++ b/org/qmckl_sherman_morrison_woodbury.org @@ -27,7 +27,7 @@ int main() { #+end_src * Naïve Sherman-Morrison - + ** ~qmckl_sherman_morrison~ :PROPERTIES: :Name: qmckl_sherman_morrison @@ -71,7 +71,7 @@ int main() { | 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 | - | double* | determinant | inout | det(Slater) the determinant of the Slater-matrix | + | double* | determinant | inout | Determinant of the Slater-matrix | *** Requirements @@ -252,7 +252,7 @@ double det = 3.407025646103221e-10; rc = qmckl_sherman_morrison(context, Dim, N_updates1, Updates1, Updates_index1, breakdown, Slater_inv1, &det); // Check that the determinant is updated properly -assert(fabs(det - 4.120398385068217e-10) < 1e-8); +assert(fabs(det + 4.120398385068217e-10) < 1e-15); for (unsigned int i = 0; i < Dim; i++) { for (unsigned int j = 0; j < Dim; j++) { @@ -298,6 +298,10 @@ assert(rc == QMCKL_SUCCESS); $B := 1 + VC$, the $2 \times 2$ matrix that is going to be inverted $D := VS^{-1}$, a $2 \times Dim$ matrix + 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_woodbury_2_args | qmckl_context | context | in | Global state | @@ -306,6 +310,7 @@ assert(rc == QMCKL_SUCCESS); | 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 | + | double* | determinant | inout | Determinant of Slater-matrix | *** Requirements @@ -328,7 +333,8 @@ assert(rc == QMCKL_SUCCESS); const double* Updates, const uint64_t* Updates_index, const double breakdown, - double* Slater_inv); + double* Slater_inv, + double* determinant); #+end_src *** C source @@ -343,7 +349,8 @@ qmckl_exit_code qmckl_woodbury_2(const qmckl_context context, const double* Updates, const uint64_t* Updates_index, const double breakdown, - double* Slater_inv) { + double* Slater_inv, + double* determinant) { /* C := S^{-1} * U, dim x 2 B := 1 + V * C, 2 x 2 @@ -381,6 +388,10 @@ qmckl_exit_code qmckl_woodbury_2(const qmckl_context context, return QMCKL_FAILURE; } + // Update det(S) when passed + if (determinant != NULL) + *determinant *= det; + // Compute B^{-1} with explicit formula for 2x2 inversion double Binv[4], idet = 1.0 / det; Binv[0] = idet * B3; @@ -427,7 +438,7 @@ qmckl_exit_code qmckl_woodbury_2(const qmckl_context context, #+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) & + (context, Dim, Updates, Updates_index, breakdown, Slater_inv, determinant) & bind(C) use, intrinsic :: iso_c_binding import @@ -439,6 +450,7 @@ qmckl_exit_code qmckl_woodbury_2(const qmckl_context context, integer (c_int64_t) , intent(in) :: Updates_index(2) real (c_double ) , intent(in) , value :: breakdown real (c_double ) , intent(inout) :: Slater_inv(Dim*Dim) + real (c_double ) , intent(inout) :: determinant end function qmckl_woodbury_2 end interface @@ -450,7 +462,9 @@ qmckl_exit_code qmckl_woodbury_2(const qmckl_context context, assert(Updates2 != NULL); assert(Updates_index2 != NULL); assert(Slater_inv2 != NULL); -rc = qmckl_woodbury_2(context, Dim, Updates2, Updates_index2, breakdown, Slater_inv2); +det = -1.4432116661319376e-11; +rc = qmckl_woodbury_2(context, Dim, Updates2, Updates_index2, breakdown, Slater_inv2, &det); +assert(fabs(det-2.367058141251457e-10) < 1e-15); for (unsigned int i = 0; i < Dim; i++) { for (unsigned int j = 0; j < Dim; j++) { res[i * Dim + j] = 0; @@ -930,7 +944,7 @@ qmckl_exit_code qmckl_sherman_morrison_smw32s(const qmckl_context context, if (remainder == 2) { 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(context, Dim, Updates_2block, Updates_index_2block, breakdown, Slater_inv); + rc = qmckl_woodbury_2(context, Dim, Updates_2block, Updates_index_2block, breakdown, Slater_inv, NULL); if (rc != 0) { // Send the entire block to slagel_splitting uint64_t l = 0; (void) qmckl_slagel_splitting(Dim, 2, Updates_2block, Updates_index_2block,