TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-10-20 14:58:23 +02:00
Code Issues Projects Releases Wiki Activity

Update determinant in Woodbury 2x2 and fix tests

Browse Source
This commit is contained in:
Pablo Oliveira 2021-10-13 10:45:54 +02:00
parent b05390a273
commit 94e9b13963
1 changed files with 22 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
org/qmckl_sherman_morrison_woodbury.org
View File

@ -27,7 +27,7 @@ int main() {
#+end_src #+end_src
* Naïve Sherman-Morrison * Naïve Sherman-Morrison
** ~qmckl_sherman_morrison~ ** ~qmckl_sherman_morrison~
:PROPERTIES: :PROPERTIES:
:Name: qmckl_sherman_morrison :Name: qmckl_sherman_morrison
@ -71,7 +71,7 @@ int main() {
| uint64_t | Updates_index[N_updates] | in | Array containing the rank-1 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 | 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 | 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 *** 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); rc = qmckl_sherman_morrison(context, Dim, N_updates1, Updates1, Updates_index1, breakdown, Slater_inv1, &det);
// Check that the determinant is updated properly // 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 i = 0; i < Dim; i++) {
for (unsigned int j = 0; j < Dim; j++) { 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 $B := 1 + VC$, the $2 \times 2$ matrix that is going to be inverted
$D := VS^{-1}$, a $2 \times Dim$ matrix $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 #+NAME: qmckl_woodbury_2_args
| qmckl_context | context | in | Global state | | 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 | | 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 | 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 | Slater_inv[Dim*Dim] | inout | Array containing the inverse of a Slater-matrix |
| double* | determinant | inout | Determinant of Slater-matrix |
*** Requirements *** Requirements
@ -328,7 +333,8 @@ assert(rc == QMCKL_SUCCESS);
const double* Updates, const double* Updates,
const uint64_t* Updates_index, const uint64_t* Updates_index,
const double breakdown, const double breakdown,
double* Slater_inv); double* Slater_inv,
double* determinant);
#+end_src #+end_src
*** C source *** C source
@ -343,7 +349,8 @@ qmckl_exit_code qmckl_woodbury_2(const qmckl_context context,
const double* Updates, const double* Updates,
const uint64_t* Updates_index, const uint64_t* Updates_index,
const double breakdown, const double breakdown,
double* Slater_inv) { double* Slater_inv,
double* determinant) {
/* /*
C := S^{-1} * U, dim x 2 C := S^{-1} * U, dim x 2
B := 1 + V * C, 2 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; return QMCKL_FAILURE;
} }
// Update det(S) when passed
if (determinant != NULL)
*determinant *= det;
// Compute B^{-1} with explicit formula for 2x2 inversion // Compute B^{-1} with explicit formula for 2x2 inversion
double Binv[4], idet = 1.0 / det; double Binv[4], idet = 1.0 / det;
Binv[0] = idet * B3; 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 #+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface interface
integer(c_int32_t) function qmckl_woodbury_2 & 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) bind(C)
use, intrinsic :: iso_c_binding use, intrinsic :: iso_c_binding
import import
@ -439,6 +450,7 @@ qmckl_exit_code qmckl_woodbury_2(const qmckl_context context,
integer (c_int64_t) , intent(in) :: Updates_index(2) integer (c_int64_t) , intent(in) :: Updates_index(2)
real (c_double ) , intent(in) , value :: breakdown real (c_double ) , intent(in) , value :: breakdown
real (c_double ) , intent(inout) :: Slater_inv(Dim*Dim) real (c_double ) , intent(inout) :: Slater_inv(Dim*Dim)
real (c_double ) , intent(inout) :: determinant
end function qmckl_woodbury_2 end function qmckl_woodbury_2
end interface end interface
@ -450,7 +462,9 @@ qmckl_exit_code qmckl_woodbury_2(const qmckl_context context,
assert(Updates2 != NULL); assert(Updates2 != NULL);
assert(Updates_index2 != NULL); assert(Updates_index2 != NULL);
assert(Slater_inv2 != 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 i = 0; i < Dim; i++) {
for (unsigned int j = 0; j < Dim; j++) { for (unsigned int j = 0; j < Dim; j++) {
res[i * Dim + j] = 0; res[i * Dim + j] = 0;
@ -930,7 +944,7 @@ qmckl_exit_code qmckl_sherman_morrison_smw32s(const qmckl_context context,
if (remainder == 2) { if (remainder == 2) {
const double *Updates_2block = &Updates[n_of_3blocks * length_3block]; const double *Updates_2block = &Updates[n_of_3blocks * length_3block];
const uint64_t *Updates_index_2block = &Updates_index[3 * n_of_3blocks]; 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 if (rc != 0) { // Send the entire block to slagel_splitting
uint64_t l = 0; uint64_t l = 0;
(void) qmckl_slagel_splitting(Dim, 2, Updates_2block, Updates_index_2block, (void) qmckl_slagel_splitting(Dim, 2, Updates_2block, Updates_index_2block,