TREX/qmckl
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-09-27 20:11:51 +02:00
Code Issues Projects Releases Wiki Activity

Merge branch 'master' into trexio

Browse Source
This commit is contained in:
Anthony Scemama 2021-10-06 23:44:19 +02:00
commit e5feaf9d0d
1 changed files with 62 additions and 42 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

104
org/qmckl_sherman_morrison_woodbury.org
View File

@ -26,9 +26,9 @@ int main() {
qmckl_exit_code rc; qmckl_exit_code rc;
#+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() {
| 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 |
*** Requirements *** Requirements
* ~context~ is not ~QMCKL_NULL_CONTEXT~ * ~context~ is not ~QMCKL_NULL_CONTEXT~
* ~Dim >= 2~ * ~Dim >= 2~
* ~N_updates >= 1~ * ~N_updates >= 1~
@ -79,7 +79,7 @@ int main() {
* ~Updates_index~ is allocated with $N_updates$ elements * ~Updates_index~ is allocated with $N_updates$ elements
* ~breakdown~ is a small number such that $0 < breakdown << 1$ * ~breakdown~ is a small number such that $0 < breakdown << 1$
* ~Slater_inv~ is allocated with $Dim \times Dim$ elements * ~Slater_inv~ is allocated with $Dim \times Dim$ elements
*** C header *** C header
#+CALL: generate_c_header(table=qmckl_sherman_morrison_args,rettyp=get_value("CRetType"),fname=get_value("Name")) #+CALL: generate_c_header(table=qmckl_sherman_morrison_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
@ -93,7 +93,7 @@ int main() {
const double* Updates, const double* Updates,
const uint64_t* Updates_index, const uint64_t* Updates_index,
const double* breakdown_p, const double* breakdown_p,
double* Slater_inv); double* Slater_inv);
#+end_src #+end_src
*** Source *** Source
@ -103,14 +103,18 @@ int main() {
#include <math.h> #include <math.h>
#include "qmckl.h" #include "qmckl.h"
qmckl_exit_code qmckl_sherman_morrison_c_(const qmckl_context context, qmckl_exit_code qmckl_sherman_morrison_c_(const qmckl_context context,
const uint64_t* Dim_p, const uint64_t* Dim_p,
const uint64_t* N_updates_p, const uint64_t* N_updates_p,
const double* Updates, const double* Updates,
const uint64_t* Updates_index, const uint64_t* Updates_index,
const double* breakdown_p, const double* breakdown_p,
double* Slater_inv) { double* Slater_inv) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
const uint64_t Dim = *Dim_p; const uint64_t Dim = *Dim_p;
const uint64_t N_updates = *N_updates_p; const uint64_t N_updates = *N_updates_p;
const double breakdown = *breakdown_p; const double breakdown = *breakdown_p;
@ -168,7 +172,7 @@ qmckl_exit_code qmckl_sherman_morrison_c_(const qmckl_context context,
:CRetType: qmckl_exit_code :CRetType: qmckl_exit_code
:FRetType: qmckl_exit_code :FRetType: qmckl_exit_code
:END: :END:
#+CALL: generate_c_interface(table=qmckl_sherman_morrison_args,rettyp=get_value("FRetType"),fname=get_value("Name")) #+CALL: generate_c_interface(table=qmckl_sherman_morrison_args,rettyp=get_value("FRetType"),fname=get_value("Name"))
#+RESULTS: #+RESULTS:
@ -194,7 +198,7 @@ qmckl_exit_code qmckl_sherman_morrison_c_(const qmckl_context context,
end function qmckl_sherman_morrison end function qmckl_sherman_morrison
#+end_src #+end_src
#+CALL: generate_f_interface(table=qmckl_sherman_morrison_args,rettyp=get_value("FRetType"),fname=get_value("Name")) #+CALL: generate_f_interface(table=qmckl_sherman_morrison_args,rettyp=get_value("FRetType"),fname=get_value("Name"))
#+RESULTS: #+RESULTS:
@ -311,7 +315,7 @@ assert(rc == QMCKL_SUCCESS);
$C:= S^{-1}U$, a Dim $\times 2$ matrix $C:= S^{-1}U$, a Dim $\times 2$ matrix
$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
#+NAME: qmckl_woodbury_2_args #+NAME: qmckl_woodbury_2_args
| qmckl_context | context | in | Global state | | qmckl_context | context | in | Global state |
@ -342,7 +346,7 @@ assert(rc == QMCKL_SUCCESS);
const double* Updates, const double* Updates,
const uint64_t* Updates_index, const uint64_t* Updates_index,
const double* breakdown_p, const double* breakdown_p,
double* Slater_inv); double* Slater_inv);
#+end_src #+end_src
*** Source *** Source
@ -352,7 +356,7 @@ assert(rc == QMCKL_SUCCESS);
#include <math.h> #include <math.h>
#include "qmckl.h" #include "qmckl.h"
qmckl_exit_code qmckl_woodbury_2_c_(const qmckl_context context, qmckl_exit_code qmckl_woodbury_2_c_(const qmckl_context context,
const uint64_t* Dim_p, const uint64_t* Dim_p,
const double* Updates, const double* Updates,
const uint64_t* Updates_index, const uint64_t* Updates_index,
@ -364,6 +368,10 @@ qmckl_exit_code qmckl_woodbury_2_c_(const qmckl_context context,
D := V * S^{-1}, 2 x dim D := V * S^{-1}, 2 x dim
,*/ ,*/
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
const uint64_t Dim = *Dim_p; const uint64_t Dim = *Dim_p;
const double breakdown = *breakdown_p; const double breakdown = *breakdown_p;
@ -514,7 +522,7 @@ assert(rc == QMCKL_SUCCESS);
#+end_src #+end_src
* Woodbury 3x3 * Woodbury 3x3
** ~qmckl_woodbury_3~ ** ~qmckl_woodbury_3~
:PROPERTIES: :PROPERTIES:
:Name: qmckl_woodbury_3 :Name: qmckl_woodbury_3
@ -559,7 +567,7 @@ assert(rc == QMCKL_SUCCESS);
const double* Updates, const double* Updates,
const uint64_t* Updates_index, const uint64_t* Updates_index,
const double* breakdown_p, const double* breakdown_p,
double* Slater_inv); double* Slater_inv);
#+end_src #+end_src
*** Source *** Source
@ -569,7 +577,7 @@ assert(rc == QMCKL_SUCCESS);
#include <math.h> #include <math.h>
#include "qmckl.h" #include "qmckl.h"
qmckl_exit_code qmckl_woodbury_3_c_(const qmckl_context context, qmckl_exit_code qmckl_woodbury_3_c_(const qmckl_context context,
const uint64_t* Dim_p, const uint64_t* Dim_p,
const double* Updates, const double* Updates,
const uint64_t* Updates_index, const uint64_t* Updates_index,
@ -581,6 +589,10 @@ qmckl_exit_code qmckl_woodbury_3_c_(const qmckl_context context,
D := V * S^{-1}, 3 x dim D := V * S^{-1}, 3 x dim
,*/ ,*/
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
const uint64_t Dim = *Dim_p; const uint64_t Dim = *Dim_p;
const double breakdown = *breakdown_p; const double breakdown = *breakdown_p;
@ -655,7 +667,7 @@ qmckl_exit_code qmckl_woodbury_3_c_(const qmckl_context context,
#+end_src #+end_src
*** Performance... *** Performance...
This function is most efficient when used in cases where there are only 3 rank-1 updates and This function is most efficient when used in cases where there are only 3 rank-1 updates and
it is sure they will not result in a singular matrix. it is sure they will not result in a singular matrix.
@ -746,7 +758,7 @@ assert(rc == QMCKL_SUCCESS);
#+end_src #+end_src
* Sherman-Morrison with update splitting * Sherman-Morrison with update splitting
** ~qmckl_sherman_morrison_splitting~ ** ~qmckl_sherman_morrison_splitting~
:PROPERTIES: :PROPERTIES:
:Name: qmckl_sherman_morrison_splitting :Name: qmckl_sherman_morrison_splitting
@ -783,7 +795,7 @@ assert(rc == QMCKL_SUCCESS);
* ~Updates_index~ is allocated with $N_updates$ elements * ~Updates_index~ is allocated with $N_updates$ elements
* ~breakdown~ is a small number such that $0 < breakdown << 1$ * ~breakdown~ is a small number such that $0 < breakdown << 1$
* ~Slater_inv~ is allocated with $Dim \times Dim$ elements * ~Slater_inv~ is allocated with $Dim \times Dim$ elements
*** C header *** C header
#+CALL: generate_c_header(table=qmckl_sherman_morrison_splitting_args,rettyp=get_value("CRetType"),fname=get_value("Name")) #+CALL: generate_c_header(table=qmckl_sherman_morrison_splitting_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
@ -797,7 +809,7 @@ 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);
#+end_src #+end_src
*** Source *** Source
@ -806,23 +818,27 @@ assert(rc == QMCKL_SUCCESS);
#include <stdbool.h> #include <stdbool.h>
#include "qmckl.h" #include "qmckl.h"
qmckl_exit_code qmckl_sherman_morrison_splitting_c_(const qmckl_context context, qmckl_exit_code qmckl_sherman_morrison_splitting_c_(const qmckl_context context,
const uint64_t* Dim, const uint64_t* Dim,
const uint64_t* N_updates, const uint64_t* N_updates,
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) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
double later_updates[*Dim * *N_updates]; double later_updates[*Dim * *N_updates];
uint64_t later_index[*N_updates]; uint64_t later_index[*N_updates];
uint64_t later = 0; uint64_t later = 0;
(void) qmckl_slagel_splitting_c(*Dim, *N_updates, Updates, Updates_index, (void) qmckl_slagel_splitting_c(*Dim, *N_updates, Updates, Updates_index,
*breakdown, Slater_inv, later_updates, later_index, &later); *breakdown, Slater_inv, later_updates, later_index, &later);
if (later > 0) { if (later > 0) {
(void) qmckl_sherman_morrison_splitting_c_(context, Dim, &later, (void) qmckl_sherman_morrison_splitting_c_(context, Dim, &later,
later_updates, later_index, breakdown, Slater_inv); later_updates, later_index, breakdown, Slater_inv);
} }
@ -834,7 +850,7 @@ qmckl_exit_code qmckl_sherman_morrison_splitting_c_(const qmckl_context context,
*** Performance... *** Performance...
This kernel performs best when there are 2 or more rank-1 update cycles and fail-rate is high. This kernel performs best when there are 2 or more rank-1 update cycles and fail-rate is high.
** C interface :noexport: ** C interface :noexport:
:PROPERTIES: :PROPERTIES:
:Name: qmckl_sherman_morrison_splitting :Name: qmckl_sherman_morrison_splitting
@ -926,7 +942,7 @@ assert(rc == QMCKL_SUCCESS);
#+end_src #+end_src
* Woodbury 3x3 and 2x2 with Sherman-Morrison and update splitting * Woodbury 3x3 and 2x2 with Sherman-Morrison and update splitting
** ~qmckl_sherman_morrison_smw32s~ ** ~qmckl_sherman_morrison_smw32s~
:PROPERTIES: :PROPERTIES:
:Name: qmckl_sherman_morrison_smw32s :Name: qmckl_sherman_morrison_smw32s
@ -972,7 +988,7 @@ assert(rc == QMCKL_SUCCESS);
const double* Updates, const double* Updates,
const uint64_t* Updates_index, const uint64_t* Updates_index,
const double* breakdown_p, const double* breakdown_p,
double* Slater_inv); double* Slater_inv);
#+end_src #+end_src
*** Source *** Source
@ -981,7 +997,7 @@ assert(rc == QMCKL_SUCCESS);
#include <stdbool.h> #include <stdbool.h>
#include "qmckl.h" #include "qmckl.h"
qmckl_exit_code qmckl_sherman_morrison_smw32s_c_(const qmckl_context context, qmckl_exit_code qmckl_sherman_morrison_smw32s_c_(const qmckl_context context,
const uint64_t* Dim_p, const uint64_t* Dim_p,
const uint64_t* N_updates_p, const uint64_t* N_updates_p,
const double* Updates, const double* Updates,
@ -989,6 +1005,10 @@ qmckl_exit_code qmckl_sherman_morrison_smw32s_c_(const qmckl_context context,
const double* breakdown_p, const double* breakdown_p,
double* Slater_inv) { double* Slater_inv) {
if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
return QMCKL_NULL_CONTEXT;
}
const uint64_t Dim = *Dim_p; const uint64_t Dim = *Dim_p;
const uint64_t N_updates = *N_updates_p; const uint64_t N_updates = *N_updates_p;
const double breakdown = *breakdown_p; const double breakdown = *breakdown_p;
@ -1141,10 +1161,10 @@ for (unsigned int i = 0; i < Dim; i++) {
} }
assert(rc == QMCKL_SUCCESS); assert(rc == QMCKL_SUCCESS);
#+end_src #+end_src
* Helper Functions * Helper Functions
Private helper-functions that are used by the Sherman-Morrison-Woodbury kernels. Private helper-functions that are used by the Sherman-Morrison-Woodbury kernels.
These functions can only be used internally by the kernels in this module. These functions can only be used internally by the kernels in this module.
** ~qmckl_slagel_splitting~ ** ~qmckl_slagel_splitting~
@ -1158,11 +1178,11 @@ These functions can only be used internally by the kernels in this module.
It is used internally to apply a collection of $N$ rank-1 updates to the inverse Slater-matrix $S^{-1}$ and It is used internally to apply a collection of $N$ rank-1 updates to the inverse Slater-matrix $S^{-1}$ and
splitting an update in two equal pieces if necessary. In case of a split, it applies the first half of the update, splitting an update in two equal pieces if necessary. In case of a split, it applies the first half of the update,
while putting the second half in a waiting queue to be applied at the end. while putting the second half in a waiting queue to be applied at the end.
Therefore, when $1+v_j^TS^{-1}u_j \geq \epsilon$, the update is applied as usual. Otherwise, $u_j$ will be redefined Therefore, when $1+v_j^TS^{-1}u_j \geq \epsilon$, the update is applied as usual. Otherwise, $u_j$ will be redefined
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}.
#+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 |
@ -1173,7 +1193,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 |
*** Requirements *** Requirements
@ -1185,8 +1205,8 @@ These functions can only be used internally by the kernels in this module.
- ~Slater_inv~ is allocated with $Dim \times Dim$ elements - ~Slater_inv~ is allocated with $Dim \times Dim$ elements
- ~later_updates~ is allocated with $later \times Dim$ elements - ~later_updates~ is allocated with $later \times Dim$ elements
- ~later_index~ is allocated with $N_updates$ elements - ~later_index~ is allocated with $N_updates$ elements
- ~later >= 0~ - ~later >= 0~
*** C header *** C header
#+CALL: generate_c_header(table=qmckl_slagel_splitting_args,rettyp=get_value("CRetType"),fname=get_value("Name")) #+CALL: generate_c_header(table=qmckl_slagel_splitting_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
@ -1202,7 +1222,7 @@ 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 );
#+end_src #+end_src
*** Source *** Source
@ -1214,11 +1234,11 @@ These functions can only be used internally by the kernels in this module.
qmckl_exit_code qmckl_slagel_splitting_c(uint64_t Dim, qmckl_exit_code qmckl_slagel_splitting_c(uint64_t Dim,
uint64_t N_updates, uint64_t N_updates,
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 *later_updates, double *later_updates,
uint64_t *later_index, uint64_t *later_index,
uint64_t *later) { uint64_t *later) {
// #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.
@ -1274,7 +1294,7 @@ qmckl_exit_code qmckl_slagel_splitting_c(uint64_t Dim,
} }
#+end_src #+end_src
*** Performance *** Performance
This function cannot be used by itself and is used in Sherman-Morrison with update splitting and Woodbury 3x3 and 2x2 This function cannot be used by itself and is used in Sherman-Morrison with update splitting and Woodbury 3x3 and 2x2
@ -1286,9 +1306,9 @@ with Sherman-Morrison and update splitting. Please look at the performance recco
:CRetType: double :CRetType: double
:FRetType: double precision :FRetType: double precision
:END: :END:
#+CALL: generate_c_interface(table=qmckl_slagel_splitting_args,rettyp=get_value("FRetType"),fname=get_value("Name")) #+CALL: generate_c_interface(table=qmckl_slagel_splitting_args,rettyp=get_value("FRetType"),fname=get_value("Name"))
*** Test :noexport: *** Test :noexport:
This kernel does not have an explicit test because it is only used internally by higher-level Sherman-Morrison-Woodbury kernels. This kernel does not have an explicit test because it is only used internally by higher-level Sherman-Morrison-Woodbury kernels.
@ -1298,7 +1318,7 @@ with Sherman-Morrison and update splitting. Please look at the performance recco
#+begin_src c :comments link :tangle (eval c_test) #+begin_src c :comments link :tangle (eval c_test)
assert (qmckl_context_destroy(context) == QMCKL_SUCCESS); assert (qmckl_context_destroy(context) == QMCKL_SUCCESS);
return 0; return 0;
} }
#+end_src #+end_src