From af54e7a7dca8ca3052e0d1e118dc660e54208f1f Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 6 Oct 2021 23:44:06 +0200 Subject: [PATCH] Checking context in SM --- org/qmckl_sherman_morrison_woodbury.org | 104 ++++++++++++++---------- 1 file changed, 62 insertions(+), 42 deletions(-) diff --git a/org/qmckl_sherman_morrison_woodbury.org b/org/qmckl_sherman_morrison_woodbury.org index cf8f15f..6fb9db0 100644 --- a/org/qmckl_sherman_morrison_woodbury.org +++ b/org/qmckl_sherman_morrison_woodbury.org @@ -26,9 +26,9 @@ int main() { qmckl_exit_code rc; #+end_src - + * Naïve Sherman-Morrison - + ** ~qmckl_sherman_morrison~ :PROPERTIES: :Name: qmckl_sherman_morrison @@ -71,7 +71,7 @@ int main() { | double | Slater_inv[Dim*Dim] | inout | Array containing the inverse of a Slater-matrix | *** Requirements - + * ~context~ is not ~QMCKL_NULL_CONTEXT~ * ~Dim >= 2~ * ~N_updates >= 1~ @@ -79,7 +79,7 @@ int main() { * ~Updates_index~ is allocated with $N_updates$ elements * ~breakdown~ is a small number such that $0 < breakdown << 1$ * ~Slater_inv~ is allocated with $Dim \times Dim$ elements - + *** C header #+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 uint64_t* Updates_index, const double* breakdown_p, - double* Slater_inv); + double* Slater_inv); #+end_src *** Source @@ -103,14 +103,18 @@ int main() { #include #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* N_updates_p, const double* Updates, const uint64_t* Updates_index, const double* breakdown_p, double* Slater_inv) { - + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + const uint64_t Dim = *Dim_p; const uint64_t N_updates = *N_updates_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 :FRetType: qmckl_exit_code :END: - + #+CALL: generate_c_interface(table=qmckl_sherman_morrison_args,rettyp=get_value("FRetType"),fname=get_value("Name")) #+RESULTS: @@ -194,7 +198,7 @@ qmckl_exit_code qmckl_sherman_morrison_c_(const qmckl_context context, end function qmckl_sherman_morrison #+end_src - + #+CALL: generate_f_interface(table=qmckl_sherman_morrison_args,rettyp=get_value("FRetType"),fname=get_value("Name")) #+RESULTS: @@ -311,7 +315,7 @@ assert(rc == QMCKL_SUCCESS); $C:= S^{-1}U$, a Dim $\times 2$ matrix $B := 1 + VC$, the $2 \times 2$ matrix that is going to be inverted $D := VS^{-1}$, a $2 \times Dim$ matrix - + #+NAME: qmckl_woodbury_2_args | qmckl_context | context | in | Global state | @@ -342,7 +346,7 @@ assert(rc == QMCKL_SUCCESS); const double* Updates, const uint64_t* Updates_index, const double* breakdown_p, - double* Slater_inv); + double* Slater_inv); #+end_src *** Source @@ -352,7 +356,7 @@ assert(rc == QMCKL_SUCCESS); #include #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 double* Updates, 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 ,*/ + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + const uint64_t Dim = *Dim_p; const double breakdown = *breakdown_p; @@ -514,7 +522,7 @@ assert(rc == QMCKL_SUCCESS); #+end_src * Woodbury 3x3 - + ** ~qmckl_woodbury_3~ :PROPERTIES: :Name: qmckl_woodbury_3 @@ -559,7 +567,7 @@ assert(rc == QMCKL_SUCCESS); const double* Updates, const uint64_t* Updates_index, const double* breakdown_p, - double* Slater_inv); + double* Slater_inv); #+end_src *** Source @@ -569,7 +577,7 @@ assert(rc == QMCKL_SUCCESS); #include #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 double* Updates, 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 ,*/ + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + const uint64_t Dim = *Dim_p; const double breakdown = *breakdown_p; @@ -655,7 +667,7 @@ qmckl_exit_code qmckl_woodbury_3_c_(const qmckl_context context, #+end_src *** Performance... - + 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. @@ -746,7 +758,7 @@ assert(rc == QMCKL_SUCCESS); #+end_src * Sherman-Morrison with update splitting - + ** ~qmckl_sherman_morrison_splitting~ :PROPERTIES: :Name: qmckl_sherman_morrison_splitting @@ -783,7 +795,7 @@ assert(rc == QMCKL_SUCCESS); * ~Updates_index~ is allocated with $N_updates$ elements * ~breakdown~ is a small number such that $0 < breakdown << 1$ * ~Slater_inv~ is allocated with $Dim \times Dim$ elements - + *** C header #+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 uint64_t* Updates_index, const double* breakdown, - double* Slater_inv); + double* Slater_inv); #+end_src *** Source @@ -806,23 +818,27 @@ assert(rc == QMCKL_SUCCESS); #include #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* N_updates, const double* Updates, const uint64_t* Updates_index, const double* breakdown, double* Slater_inv) { - + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + double later_updates[*Dim * *N_updates]; uint64_t later_index[*N_updates]; uint64_t later = 0; (void) qmckl_slagel_splitting_c(*Dim, *N_updates, Updates, Updates_index, *breakdown, Slater_inv, later_updates, later_index, &later); - + 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); } @@ -834,7 +850,7 @@ qmckl_exit_code qmckl_sherman_morrison_splitting_c_(const qmckl_context context, *** Performance... This kernel performs best when there are 2 or more rank-1 update cycles and fail-rate is high. - + ** C interface :noexport: :PROPERTIES: :Name: qmckl_sherman_morrison_splitting @@ -926,7 +942,7 @@ assert(rc == QMCKL_SUCCESS); #+end_src * Woodbury 3x3 and 2x2 with Sherman-Morrison and update splitting - + ** ~qmckl_sherman_morrison_smw32s~ :PROPERTIES: :Name: qmckl_sherman_morrison_smw32s @@ -972,7 +988,7 @@ assert(rc == QMCKL_SUCCESS); const double* Updates, const uint64_t* Updates_index, const double* breakdown_p, - double* Slater_inv); + double* Slater_inv); #+end_src *** Source @@ -981,7 +997,7 @@ assert(rc == QMCKL_SUCCESS); #include #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* N_updates_p, const double* Updates, @@ -989,6 +1005,10 @@ qmckl_exit_code qmckl_sherman_morrison_smw32s_c_(const qmckl_context context, const double* breakdown_p, double* Slater_inv) { + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + const uint64_t Dim = *Dim_p; const uint64_t N_updates = *N_updates_p; const double breakdown = *breakdown_p; @@ -1141,10 +1161,10 @@ for (unsigned int i = 0; i < Dim; i++) { } assert(rc == QMCKL_SUCCESS); #+end_src - + * 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. ** ~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 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. - + 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 $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 | 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 | @@ -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 | | 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 | - + *** 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 - ~later_updates~ is allocated with $later \times Dim$ elements - ~later_index~ is allocated with $N_updates$ elements - - ~later >= 0~ - + - ~later >= 0~ + *** C header #+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* later_updates, uint64_t* later_index, - uint64_t* later ); + uint64_t* later ); #+end_src *** 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, uint64_t N_updates, - const double *Updates, + const double *Updates, const uint64_t *Updates_index, const double breakdown, double *Slater_inv, - double *later_updates, + double *later_updates, uint64_t *later_index, uint64_t *later) { // #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 - + *** Performance 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 :FRetType: double precision :END: - + #+CALL: generate_c_interface(table=qmckl_slagel_splitting_args,rettyp=get_value("FRetType"),fname=get_value("Name")) - + *** Test :noexport: 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) assert (qmckl_context_destroy(context) == QMCKL_SUCCESS); return 0; - + } #+end_src