TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-08-18 03:11:44 +02:00
Code Issues Projects Releases Wiki Activity

Update determinant in the SM+WB splitting versions

Browse Source
This commit is contained in:
Pablo Oliveira 2021-10-13 11:55:20 +02:00
parent 9ca88679f9
commit 4f2e8b6d8e
1 changed files with 49 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

65
org/qmckl_sherman_morrison_woodbury.org
View File

@ -728,6 +728,9 @@ assert(rc == QMCKL_SUCCESS);
case the Slater-matrix that would have resulted from applying the updates is singular and therefore the case the Slater-matrix that would have resulted from applying the updates is singular and therefore the
kernel exits with an exit code. kernel exits with an exit code.
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_splitting_args #+NAME: qmckl_sherman_morrison_splitting_args
| qmckl_context | context | in | Global state | | qmckl_context | context | in | Global state |
| uint64_t | Dim | in | Leading dimension of Slater_inv | | uint64_t | Dim | in | Leading dimension of Slater_inv |
@ -736,6 +739,10 @@ assert(rc == QMCKL_SUCCESS);
| 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 | Determinant of the Slater-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.
*** Requirements *** Requirements
@ -760,7 +767,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
@ -775,7 +783,8 @@ qmckl_exit_code qmckl_sherman_morrison_splitting(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) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT; return QMCKL_NULL_CONTEXT;
@ -786,11 +795,11 @@ qmckl_exit_code qmckl_sherman_morrison_splitting(const qmckl_context context,
uint64_t later = 0; uint64_t later = 0;
(void) qmckl_slagel_splitting(Dim, N_updates, Updates, Updates_index, (void) qmckl_slagel_splitting(Dim, N_updates, Updates, Updates_index,
breakdown, Slater_inv, later_updates, later_index, &later); breakdown, Slater_inv, later_updates, later_index, &later, determinant);
if (later > 0) { if (later > 0) {
(void) qmckl_sherman_morrison_splitting(context, Dim, later, (void) qmckl_sherman_morrison_splitting(context, Dim, later,
later_updates, later_index, breakdown, Slater_inv); later_updates, later_index, breakdown, Slater_inv, determinant);
} }
return QMCKL_SUCCESS; return QMCKL_SUCCESS;
@ -815,7 +824,7 @@ qmckl_exit_code qmckl_sherman_morrison_splitting(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_sherman_morrison_splitting & integer(c_int32_t) function qmckl_sherman_morrison_splitting &
(context, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv) & (context, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) &
bind(C) bind(C)
use, intrinsic :: iso_c_binding use, intrinsic :: iso_c_binding
import import
@ -828,6 +837,7 @@ qmckl_exit_code qmckl_sherman_morrison_splitting(const qmckl_context context,
integer (c_int64_t) , intent(in) :: Updates_index(N_updates) integer (c_int64_t) , intent(in) :: Updates_index(N_updates)
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_sherman_morrison_splitting end function qmckl_sherman_morrison_splitting
end interface end interface
@ -839,7 +849,9 @@ qmckl_exit_code qmckl_sherman_morrison_splitting(const qmckl_context context,
assert(Updates3 != NULL); assert(Updates3 != NULL);
assert(Updates_index3 != NULL); assert(Updates_index3 != NULL);
assert(Slater_inv3_2 != NULL); assert(Slater_inv3_2 != NULL);
rc = qmckl_sherman_morrison_splitting(context, Dim, N_updates3, Updates3, Updates_index3, breakdown, Slater_inv3_2); det = -1.23743195512859e-09;
rc = qmckl_sherman_morrison_splitting(context, Dim, N_updates3, Updates3, Updates_index3, breakdown, Slater_inv3_2, &det);
assert(fabs(det - 1.602708950725074e-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;
@ -877,6 +889,9 @@ assert(rc == QMCKL_SUCCESS);
with Woodbury 2x2 instead of sending them all to Sherman-Morrison with update splitting. For example, in the case of with Woodbury 2x2 instead of sending them all to Sherman-Morrison with update splitting. For example, in the case of
5 updates the updates are applied in 1 block of 3 updates end 1 block of 2 updates. 5 updates the updates are applied in 1 block of 3 updates end 1 block of 2 updates.
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_smw32s_args #+NAME: qmckl_sherman_morrison_smw32s_args
| qmckl_context | context | in | Global state | | qmckl_context | context | in | Global state |
| uint64_t | Dim | in | Leading dimension of Slater_inv | | uint64_t | Dim | in | Leading dimension of Slater_inv |
@ -885,6 +900,8 @@ assert(rc == QMCKL_SUCCESS);
| 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 | Determinant of the Slater-matrix |
*** Requirements *** Requirements
@ -909,7 +926,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
@ -924,7 +942,8 @@ qmckl_exit_code qmckl_sherman_morrison_smw32s(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) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT; return QMCKL_NULL_CONTEXT;
@ -945,11 +964,11 @@ qmckl_exit_code qmckl_sherman_morrison_smw32s(const qmckl_context context,
for (uint64_t i = 0; i < n_of_3blocks; i++) { for (uint64_t i = 0; i < n_of_3blocks; i++) {
const double *Updates_3block = &Updates[i * length_3block]; const double *Updates_3block = &Updates[i * length_3block];
const uint64_t *Updates_index_3block = &Updates_index[i * 3]; const uint64_t *Updates_index_3block = &Updates_index[i * 3];
rc = qmckl_woodbury_3(context, Dim, Updates_3block, Updates_index_3block, breakdown, Slater_inv, NULL); rc = qmckl_woodbury_3(context, Dim, Updates_3block, Updates_index_3block, breakdown, Slater_inv, determinant);
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;
rc = qmckl_slagel_splitting(Dim, 3, Updates_3block, Updates_index_3block, rc = qmckl_slagel_splitting(Dim, 3, Updates_3block, Updates_index_3block,
breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l); breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l, determinant);
later = later + l; later = later + l;
} }
} }
@ -959,11 +978,11 @@ 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, NULL); rc = qmckl_woodbury_2(context, Dim, Updates_2block, Updates_index_2block, breakdown, Slater_inv, determinant);
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,
breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l); breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l, determinant);
later = later + l; later = later + l;
} }
} }
@ -973,12 +992,12 @@ qmckl_exit_code qmckl_sherman_morrison_smw32s(const qmckl_context context,
const uint64_t *Updates_index_1block = &Updates_index[3 * n_of_3blocks]; const uint64_t *Updates_index_1block = &Updates_index[3 * n_of_3blocks];
uint64_t l = 0; uint64_t l = 0;
(void) qmckl_slagel_splitting(Dim, 1, Updates_1block, Updates_index_1block, (void) qmckl_slagel_splitting(Dim, 1, Updates_1block, Updates_index_1block,
breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l); breakdown, Slater_inv, later_updates + (Dim * later), later_index + later, &l, determinant);
later = later + l; later = later + l;
} }
if (later > 0) { if (later > 0) {
(void) qmckl_sherman_morrison_splitting(context, Dim, later, later_updates, later_index, breakdown, Slater_inv); (void) qmckl_sherman_morrison_splitting(context, Dim, later, later_updates, later_index, breakdown, Slater_inv, determinant);
} }
return QMCKL_SUCCESS; return QMCKL_SUCCESS;
} }
@ -1026,7 +1045,10 @@ qmckl_exit_code qmckl_sherman_morrison_smw32s(const qmckl_context context,
assert(Updates5 != NULL); assert(Updates5 != NULL);
assert(Updates_index5 != NULL); assert(Updates_index5 != NULL);
assert(Slater_inv5 != NULL); assert(Slater_inv5 != NULL);
rc = qmckl_sherman_morrison_smw32s(context, Dim, N_updates5, Updates5, Updates_index5, breakdown, Slater_inv5); det = -3.186005284713128e-10;
rc = qmckl_sherman_morrison_smw32s(context, Dim, N_updates5, Updates5, Updates_index5, breakdown, Slater_inv5, &det);
assert(fabs(det + 5.260200118412903e-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;
@ -1070,6 +1092,9 @@ These functions can only be used internally by the kernels in this module.
as $\frac{1}{2}u_j$. One half is applied immediately, the other half will be applied at the end of the algorithm, using vectors as $\frac{1}{2}u_j$. One half is applied immediately, the other half will be applied at the end of the algorithm, using vectors
$u_{j'}=\frac{1}{2}u_j$ and $v_{j'}^T=v_{j}^T$, which are stored in the array \texttt{later_updates}. $u_{j'}=\frac{1}{2}u_j$ and $v_{j'}^T=v_{j}^T$, which are stored in the array \texttt{later_updates}.
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_slagel_splitting_args #+NAME: qmckl_slagel_splitting_args
| uint64_t | Dim | in | Leading dimension of Slater_inv | | 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 | | uint64_t | N_updates | in | Number of rank-1 updates to be applied to Slater_inv |
@ -1080,6 +1105,7 @@ These functions can only be used internally by the kernels in this module.
| double | later_updates[Dim * N_updates] | inout | Array containing the split updates for later | | double | later_updates[Dim * N_updates] | inout | Array containing the split updates for later |
| uint64_t | later_index[N_updates] | inout | Array containing the positions of the split updates for later | | uint64_t | later_index[N_updates] | inout | Array containing the positions of the split updates for later |
| uint64_t | later | inout | Number of split updates for later | | uint64_t | later | inout | Number of split updates for later |
| double* | determinant | inout | Determinant of the Slater-matrix |
*** Requirements *** Requirements
@ -1109,7 +1135,8 @@ These functions can only be used internally by the kernels in this module.
double* Slater_inv, double* Slater_inv,
double* later_updates, double* later_updates,
uint64_t* later_index, uint64_t* later_index,
uint64_t* later ); uint64_t* later,
double* determinant );
#+end_src #+end_src
*** C source *** C source
@ -1127,7 +1154,8 @@ qmckl_exit_code qmckl_slagel_splitting(uint64_t Dim,
double *Slater_inv, double *Slater_inv,
double *later_updates, double *later_updates,
uint64_t *later_index, uint64_t *later_index,
uint64_t *later) { uint64_t *later,
double *determinant) {
// #ifdef DEBUG // Leave commented out since debugging information is not yet implemented in QMCkl. // #ifdef DEBUG // Leave commented out since debugging information is not yet implemented in QMCkl.
// std::cerr << "Called slagel_splitting with " << N_updates << " updates" << std::endl; // std::cerr << "Called slagel_splitting with " << N_updates << " updates" << std::endl;
// #endif // #endif
@ -1162,6 +1190,9 @@ qmckl_exit_code qmckl_slagel_splitting(uint64_t Dim,
} // From here onwards we continue with applying the first havel of the update to Slater_inv } // From here onwards we continue with applying the first havel of the update to Slater_inv
double iden = 1 / den; double iden = 1 / den;
if (determinant != NULL)
*determinant *= den;
// D = v^T x S^{-1} // D = v^T x S^{-1}
for (uint64_t j = 0; j < Dim; j++) { for (uint64_t j = 0; j < Dim; j++) {
D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j]; D[j] = Slater_inv[(Updates_index[l] - 1) * Dim + j];