TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2025-01-03 18:16:28 +01:00
Code Issues Projects Releases Wiki Activity

Update determinant in Woodbury 3x3

Browse Source
This commit is contained in:
Pablo Oliveira 2021-10-13 10:51:42 +02:00
parent 94e9b13963
commit 9ca88679f9
1 changed files with 20 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
org/qmckl_sherman_morrison_woodbury.org
View File

@ -504,6 +504,11 @@ assert(rc == QMCKL_SUCCESS);
$B := 1 + VC$, the $3 \times 3$ matrix that is going to be inverted
$D := VS^{-1}$, a $3 \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_3_args
| qmckl_context | context | in | Global state |
| uint64_t | Dim | in | Leading dimension of Slater_inv |
@ -511,6 +516,7 @@ assert(rc == QMCKL_SUCCESS);
| uint64_t | Updates_index[3] | 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
@ -533,7 +539,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
@ -548,7 +555,8 @@ qmckl_exit_code qmckl_woodbury_3(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 3
B := 1 + V * C, 3 x 3
@ -594,6 +602,10 @@ qmckl_exit_code qmckl_woodbury_3(const qmckl_context context,
return QMCKL_FAILURE;
}
// Update det(Slater) if passed
if (determinant != NULL)
*determinant *= det;
// Compute B^{-1} with explicit formula for 3x3 inversion
double Binv[9], idet = 1.0 / det;
Binv[0] = (B4 * B8 - B7 * B5) * idet;
@ -647,7 +659,7 @@ qmckl_exit_code qmckl_woodbury_3(const qmckl_context context,
#+begin_src f90 :tangle (eval fh_func) :comments org :exports none
interface
integer(c_int32_t) function qmckl_woodbury_3 &
(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
@ -659,6 +671,7 @@ qmckl_exit_code qmckl_woodbury_3(const qmckl_context context,
integer (c_int64_t) , intent(in) :: Updates_index(3)
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_3
end interface
@ -670,7 +683,9 @@ qmckl_exit_code qmckl_woodbury_3(const qmckl_context context,
assert(Updates3 != NULL);
assert(Updates_index3 != NULL);
assert(Slater_inv3_1 != NULL);
rc = qmckl_woodbury_3(context, Dim, Updates3, Updates_index3, breakdown, Slater_inv3_1);
det = -1.23743195512859e-09;
rc = qmckl_woodbury_3(context, Dim, Updates3, Updates_index3, breakdown, Slater_inv3_1, &det);
assert(fabs(det - 1.602708950725074e-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 +945,7 @@ qmckl_exit_code qmckl_sherman_morrison_smw32s(const qmckl_context context,
for (uint64_t i = 0; i < n_of_3blocks; i++) {
const double *Updates_3block = &Updates[i * length_3block];
const uint64_t *Updates_index_3block = &Updates_index[i * 3];
rc = qmckl_woodbury_3(context, Dim, Updates_3block, Updates_index_3block, breakdown, Slater_inv);
rc = qmckl_woodbury_3(context, Dim, Updates_3block, Updates_index_3block, breakdown, Slater_inv, NULL);
if (rc != 0) { // Send the entire block to slagel_splitting
uint64_t l = 0;
rc = qmckl_slagel_splitting(Dim, 3, Updates_3block, Updates_index_3block,