From 5e5c15a09d7ebe8dd77264f6cc24aecce38e05b5 Mon Sep 17 00:00:00 2001 From: Francois Coppens Date: Mon, 20 Feb 2023 18:31:39 +0100 Subject: [PATCH 1/8] - Added qmckl_context to Slagel Splitting kernel - Renamed it to Sherman-Morrison Splitting Core. - Sherman-Morrison Splitting Core now callable on its own. User is responsible for what to do with the output data. - Added default switch cases with asserts to generate crash with message if a template for a specific size is missing. - Added switch breaks to prevent the default case to always execute and make the kernel crash at the assert. - Reorganisded the Sherman-Morrison Splitting kernel so that the HPC variant always calls the Core HPC variant and not the generec variant and make duplicate decisions. --- org/qmckl_sherman_morrison_woodbury.org | 1305 ++++++++++++----------- 1 file changed, 709 insertions(+), 596 deletions(-) diff --git a/org/qmckl_sherman_morrison_woodbury.org b/org/qmckl_sherman_morrison_woodbury.org index 98193cc..4c69d63 100644 --- a/org/qmckl_sherman_morrison_woodbury.org +++ b/org/qmckl_sherman_morrison_woodbury.org @@ -29,7 +29,7 @@ int main() { This is the range that determines the how many high performance kernel instantces will be generated, using the C-function templates defined in the sections below. If the name of the C-function template is called ~qmckl_kernel_{Dim}~, then ~range(K, L+1)~ will results in kernel instances from ~qmckl_kernel_K~ to ~qmckl_kernel_L~. #+NAME:kernel_generator_range #+begin_src python :noweb yes :exports none -range(2, 22) +range(2, 3) #+end_src @@ -272,7 +272,7 @@ end function qmckl_sherman_morrison_naive_doc #+CALL: generate_c_header(table=qmckl_sherman_morrison_naive_args,rettyp=get_value("CRetType"),fname=get_value("Name")) #+RESULTS: -#+begin_src c :tangle (eval h_func) :comments org +#+begin_src c :tangle (eval h_func) : comments org qmckl_exit_code qmckl_sherman_morrison_naive ( const qmckl_context context, const uint64_t LDS, @@ -288,7 +288,7 @@ qmckl_exit_code qmckl_sherman_morrison_naive ( #+CALL: generate_c_header(table=qmckl_sherman_morrison_naive_args,rettyp=get_value("CRetType"),fname="qmckl_sherman_morrison_naive_hpc") #+RESULTS: -#+begin_src c :tangle (eval h_func) :comments org +#+begin_src c :tangle (eval h_private_func) :comments org qmckl_exit_code qmckl_sherman_morrison_naive_hpc ( const qmckl_context context, const uint64_t LDS, @@ -324,6 +324,7 @@ Common includes and macros used by all the Sherman-Morrison-Woodbury kernels. #include #include "qmckl.h" #include "config.h" +#include "assert.h" // Order important because // __GNUC__ also set in ICC, ICX and CLANG @@ -736,7 +737,603 @@ assert(rc == QMCKL_SUCCESS); #+end_src -* Sherman-Morrison with update splitting +* Sherman-Morrison with Slagel Splitting (core) +** ~qmckl_slagel_splitting~ +:PROPERTIES: +:Name: qmckl_slagel_splitting +:CRetType: qmckl_exit_code +:FRetType: qmckl_exit_code +:END: + +*** Introduction +~qmckl_slagel_splitting~ is the inner core part of 'Sherman-Morrison with update splitting' in the next section. +It is not normally used by itself but it is possible to use it nonetheless. + +It has three extra parameters in its API: +- ~later_updates~ initially empty array that will contain the second halves of updates that were split during kernel execution +- ~later_index~ initially empty array that will contain the row/column numbers of the updates that were split during execution +- ~later~ initially zero integer that records the number of updates that were split during exection. + +It is up to the user to decide what to do with these updates once the kernel returns. Normally ~qmckl_slagel_splitting~ is +used as the core part of a recursive function, as is done in ~qmckl_sherman_morrison_splitting~ or as part of a more complex +kernel like ~qmckl_sherman_morrison_smw32s~. + +If the determinant is passed it will only be partially updated if there were any update splits. + +*** API +#+NAME: qmckl_slagel_splitting_args +| Variable | Type | In/Out | Description | +|-----------------+-------------------------+--------+---------------------------------------------------------------| +| ~context~ | ~qmckl_context~ | in | Global state | +| ~LDS~ | ~uint64_t~ | in | Leading dimension of Slater_inv | +| ~Dim~ | ~uint64_t~ | in | Dimension of Slater_inv | +| ~N_updates~ | ~uint64_t~ | in | Number of rank-1 updates to be applied to Slater_inv | +| ~Updates~ | ~double[N_updates*LDS]~ | in | Array containing the rank-1 updates | +| ~Updates_index~ | ~uint64_t[N_updates]~ | in | Array containing positions of the rank-1 updates | +| ~breakdown~ | ~double~ | in | Break-down parameter on which to fail or not | +| ~Slater_inv~ | ~double[Dim*LDS]~ | inout | Array containing the inverse Slater-matrix | +| ~later_updates~ | ~double[N_updates*LDS]~ | inout | Array containing the split updates for later | +| ~later_index~ | ~uint64_t[N_updates]~ | inout | Array containing the positions of the split updates for later | +| ~later~ | ~uint64_t~ | inout | Number of split updates for later | +| ~determinant~ | ~double~ | inout | Determinant of the Slater-matrix | + +*** Requirements + * ~LDS >= 2~ + * ~Dim >= 2~ + * ~N_updates >= 1~ + * ~Updates~ is allocated with $N_updates \times Dim$ elements + * ~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 + * ~later_updates~ is allocated with $later \times Dim$ elements + * ~later_index~ is allocated with $N_updates$ elements + * ~later >= 0~ + +*** Pedagogical kernel source (in Fortran) +The following source code written in Fortran is inteded to illustrate how the kernel works. Even though the kernel is +able to do numerically correct computations, it does not do it in the most efficient way possible. It should therefore +not be used in real workloads. + +#+begin_src f90 :tangle (eval f) +integer function qmckl_slagel_splitting_doc_f( & + context, & + lds, dim, & + nupdates, & + upds, & + updates_index, & + breakdown, & + s_inv, & + later_updates, & + later_index, & + later, & + determinant) result(info) + + use qmckl + implicit none + integer*8 , intent(in) :: context + integer*8 , intent(in) :: lds, dim + integer*8 , intent(in) :: nupdates + integer*8 , intent(in) :: updates_index(nupdates) + real*8 , intent(in) :: upds(nupdates * lds) + real*8 , intent(in) :: breakdown + real*8 , intent(inout) :: s_inv(dim * lds) + real*8 , intent(inout) :: later_updates(nupdates * lds) + integer*8 , intent(inout) :: later_index(nupdates) + integer*8 , intent(inout) :: later + real*8 , intent(inout) :: determinant + + real*8 , dimension(lds, nupdates) :: Updates + real*8 , dimension(dim, lds) :: Inverse + + info = QMCKL_FAILURE + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + ! Convert 'upds' and 's_inv' into the more easily readable Fortran + ! matrices 'Updates' and 'Inverse'. + call convert(upds, s_inv, Updates, Inverse, nupdates, lds, dim) + + ! YET TO BE IMPLEMENTED + + ! Copy updated inverse back to s_inv + call copy_back(Inverse, s_inv, lds, dim) + + info = QMCKL_SUCCESS + +end function qmckl_slagel_splitting_doc_f +#+end_src + +**** C interface to the pedagogical kernel (not directly exposed) +The following Fortran function ~qmckl_slagel_splitting_doc~ makes sure +that the pedagogical kernel ~qmckl_slagel_splitting_doc_f~, written in +Fortran, can be called from C using the ~ISO_C_BINDING~. The Fortran function +~qmckl_slagel_splitting_doc~ will be exposed in the header file 'qmckl.h' +for C users and in the module file 'qmckl_f.F90' for Fortran users. + +#+CALL: generate_c_interface(table=qmckl_slagel_splitting_args,rettyp=get_value("CRetType"),fname="qmckl_slagel_splitting_doc") + +#+RESULTS: +#+begin_src f90 :tangle (eval f) :comments org :exports none +integer(c_int32_t) function qmckl_slagel_splitting_doc & + (context, & + LDS, & + Dim, & + N_updates, & + Updates, & + Updates_index, & + breakdown, & + Slater_inv, & + later_updates, & + later_index, & + later, & + determinant) & + bind(C) result(info) + + use, intrinsic :: iso_c_binding + implicit none + + integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: LDS + integer (c_int64_t) , intent(in) , value :: Dim + integer (c_int64_t) , intent(in) , value :: N_updates + real (c_double ) , intent(in) :: Updates(N_updates*LDS) + integer (c_int64_t) , intent(in) :: Updates_index(N_updates) + real (c_double ) , intent(in) , value :: breakdown + real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS) + real (c_double ) , intent(inout) :: later_updates(N_updates*LDS) + integer (c_int64_t) , intent(inout) :: later_index(N_updates) + integer (c_int64_t) , intent(inout) :: later + real (c_double ) , intent(inout) :: determinant + + integer(c_int32_t), external :: qmckl_slagel_splitting_doc_f + info = qmckl_slagel_splitting_doc_f & + (context, & + LDS, & + Dim, & + N_updates, & + Updates, & + Updates_index, & + breakdown, & + Slater_inv, & + later_updates, & + later_index, & + later, & + determinant) + +end function qmckl_slagel_splitting_doc +#+end_src + +*** C headers (exposed in qmckl.h) +#+CALL: generate_c_header(table=qmckl_slagel_splitting_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + +#+RESULTS: +#+begin_src c :tangle (eval h_func) :comments org +qmckl_exit_code qmckl_slagel_splitting ( + const qmckl_context context, + const uint64_t LDS, + const uint64_t Dim, + const uint64_t N_updates, + const double* Updates, + const uint64_t* Updates_index, + const double breakdown, + double* Slater_inv, + double* later_updates, + uint64_t* later_index, + uint64_t* later, + double* determinant ); +#+end_src + +#+CALL: generate_c_header(table=qmckl_slagel_splitting_args,rettyp=get_value("CRetType"),fname="qmckl_slagel_splitting_doc") + +#+RESULTS: +#+begin_src c :tangle (eval h_func) :comments org +qmckl_exit_code qmckl_slagel_splitting_doc ( + const qmckl_context context, + const uint64_t LDS, + const uint64_t Dim, + const uint64_t N_updates, + const double* Updates, + const uint64_t* Updates_index, + const double breakdown, + double* Slater_inv, + double* later_updates, + uint64_t* later_index, + uint64_t* later, + double* determinant ); +#+end_src + +*** C sources +#+begin_src c :tangle (eval c) :comments org +qmckl_exit_code qmckl_slagel_splitting_hpc( + const qmckl_context context, + uint64_t LDS, + uint64_t Dim, + uint64_t N_updates, + const double* __restrict Updates, + const uint64_t* __restrict Updates_index, + const double breakdown, + double* __restrict Slater_inv, + double* __restrict later_updates, + uint64_t* __restrict later_index, + uint64_t* __restrict later, + double* __restrict determinant) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return qmckl_failwith( + context, + QMCKL_NULL_CONTEXT, + "qmckl_slagel_splitting_hpc", + NULL); + } + + double __attribute__((aligned(8))) C[LDS]; + double __attribute__((aligned(8))) D[LDS]; + + uint64_t l = 0; + // For each update + while (l < N_updates) { + // C = S^{-1} x U_l + for (uint64_t i = 0; i < Dim; i++) { + C[i] = 0.0f; + IVDEP + ALIGNED + for (uint64_t j = 0; j < LDS; j++) { + C[i] += Slater_inv[i * LDS + j] * Updates[l * LDS + j]; + } + } + + // Denominator + const int cui = Updates_index[l] - 1; + double den = 1.0f + C[cui]; + if (fabs(den) < breakdown) { + // U_l = U_l / 2: split the update in 2 equal halves and save the + // second halve in later_updates + IVDEP + ALIGNED + for (uint64_t i = 0; i < LDS; i++) { + later_updates[*later * LDS + i] = Updates[l * LDS + i] * 0.5f; + C[i] *= 0.5f; + } + later_index[*later] = Updates_index[l]; + (*later)++; + + den = 1.0f + C[cui]; + } // From here onwards we continue with applying the first halve of the + // update to Slater_inv + double iden = 1.0f / den; + + if (determinant) + *determinant *= den; + + // D = v^T x S^{-1} : 1 x LDS + IVDEP + ALIGNED + for (uint64_t j = 0; j < LDS; j++) { + D[j] = Slater_inv[cui * LDS + j]; + } + + // S^{-1} = S^{-1} - C x D / den + for (uint64_t i = 0; i < Dim; i++) { + IVDEP + ALIGNED + for (uint64_t j = 0; j < LDS; j++) { + const double update = C[i] * D[j] * iden; + Slater_inv[i * LDS + j] -= update; + } + } + l += 1; + } + + return QMCKL_SUCCESS; +} + #+end_src + +#+NAME:slagel_splitting_template_code +#+begin_src c +static inline qmckl_exit_code qmckl_slagel_splitting_{Dim}( + const qmckl_context context, + uint64_t N_updates, + const double* __restrict Updates, + const uint64_t* __restrict Updates_index, + const double breakdown, + double* __restrict Slater_inv, + double* __restrict later_updates, + uint64_t* __restrict later_index, + uint64_t* __restrict later, + double* __restrict determinant) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return qmckl_failwith( + context, + QMCKL_NULL_CONTEXT, + "qmckl_slagel_splitting_{Dim}", + NULL); + } + + double __attribute__((aligned(8))) C[D{Dim}_P]; + double __attribute__((aligned(8))) D[D{Dim}_P]; + + uint64_t l = 0; + // For each update + while (l < N_updates) { + // C = S^{-1} x U_l + for (uint64_t i = 0; i < {Dim}; i++) { + C[i] = 0.0f; + IVDEP + ALIGNED + for (uint64_t j = 0; j < D{Dim}_P; j++) { + C[i] += Slater_inv[i * D{Dim}_P + j] * Updates[l * D{Dim}_P + j]; + } + } + + // Denominator + const int cui = Updates_index[l] - 1; + double den = 1.0f + C[cui]; + if (fabs(den) < breakdown) { + // U_l = U_l / 2: split the update in 2 equal halves and save the + // second halve in later_updates + IVDEP + ALIGNED + for (uint64_t i = 0; i < D{Dim}_P; i++) { + later_updates[*later * D{Dim}_P + i] = Updates[l * D{Dim}_P + i] * 0.5f; + C[i] *= 0.5f; + } + later_index[*later] = Updates_index[l]; + (*later)++; + + den = 1.0f + C[cui]; + } // From here onwards we continue with applying the first halve of the + // update to Slater_inv + double iden = 1.0f / den; + + if (determinant) + *determinant *= den; + + // D = v^T x S^{-1} : 1 x D{Dim}_P + IVDEP + ALIGNED + for (uint64_t j = 0; j < D{Dim}_P; j++) { + D[j] = Slater_inv[cui * D{Dim}_P + j]; + } + + // S^{-1} = S^{-1} - C x D / den + for (uint64_t i = 0; i < {Dim}; i++) { + IVDEP + ALIGNED + for (uint64_t j = 0; j < D{Dim}_P; j++) { + const double update = C[i] * D[j] * iden; + Slater_inv[i * D{Dim}_P + j] -= update; + } + } + l += 1; + } + + return QMCKL_SUCCESS; +} +#+end_src + +#+NAME:slagel_splitting_kernel_generator +#+begin_src python :noweb yes +text=""" +<> +""" +result = [] +for Dim in <>: + Dim=str(Dim) + result.append(text.replace("{Dim}",Dim) ) + +return '\n'.join(result) +#+end_src + +#+NAME:slagel_splitting_switch-case_generator +#+begin_src python :noweb yes +text=""" +case {Dim}: { + return qmckl_slagel_splitting_{Dim}( + context, + N_updates, + Updates, + Updates_index, + breakdown, + Slater_inv, + later_updates, + later_index, + later, + determinant); + break; +} +""" +result = [] +for Dim in <>: + Dim=str(Dim) + result.append(text.replace("{Dim}",Dim) ) + +return '\n'.join(result) +#+end_src + +#+begin_src c :tangle (eval c) :comments org :noweb yes +<> +#+end_src + +#+begin_src c :tangle (eval c) :comments org :noweb yes +qmckl_exit_code qmckl_slagel_splitting( + const qmckl_context context, + const uint64_t LDS, + const uint64_t Dim, + const uint64_t N_updates, + const double* Updates, + const uint64_t* Updates_index, + const double breakdown, + double* Slater_inv, + double* later_updates, + uint64_t* later_index, + uint64_t* later, + double* determinant) { + + #ifdef HAVE_HPC + if (LDS == (1+(Dim-1)/SIMD_LENGTH)*SIMD_LENGTH) { // Most cases + switch (Dim) { + <> + case default: { + assert(0 == 1 && "TEMPLATE NOT IMPLEMENTED!"); + break; + } + } + } + else { // Updating smaller sub-matrix + return qmckl_slagel_splitting_hpc( + context, + LDS, + Dim, + N_updates, + Updates, + Updates_index, + breakdown, + Slater_inv, + later_updates, + later_index, + later, + determinant); + } + #else + // return qmckl_slagel_splitting_doc( + // context, + // LDS, + // Dim, + // N_updates, + // Updates, + // Updates_index, + // breakdown, + // Slater_inv, + // later_updates, + // later_index, + // later, + // determinant); + return qmckl_slagel_splitting_hpc( + context, + LDS, + Dim, + N_updates, + Updates, + Updates_index, + breakdown, + Slater_inv, + later_updates, + later_index, + later, + determinant); + #endif + + return QMCKL_FAILURE; +} +#+end_src + +*** Fortran interfaces (exposed in qmckl_f.F90) +:PROPERTIES: +:Name: qmckl_slagel_splitting +:CRetType: qmckl_exit_code +:FRetType: qmckl_exit_code +:END: + +#+CALL: generate_f_interface(table=qmckl_slagel_splitting_args,rettyp=get_value("FRetType"),fname="qmckl_slagel_splitting_hpc") + +#+RESULTS: +#+begin_src f90 :tangle (eval fh_func) :comments org :exports none +interface + integer(c_int32_t) function qmckl_slagel_splitting_hpc & + (context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, & + Slater_inv, later_updates, later_index, later, determinant) & + bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + + integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: LDS + integer (c_int64_t) , intent(in) , value :: Dim + integer (c_int64_t) , intent(in) , value :: N_updates + real (c_double ) , intent(in) :: Updates(N_updates*LDS) + integer (c_int64_t) , intent(in) :: Updates_index(N_updates) + real (c_double ) , intent(in) , value :: breakdown + real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS) + real (c_double ) , intent(inout) :: later_updates(N_updates*LDS) + integer (c_int64_t) , intent(inout) :: later_index(N_updates) + integer (c_int64_t) , intent(inout) :: later + real (c_double ) , intent(inout) :: determinant + + end function qmckl_slagel_splitting_hpc +end interface +#+end_src + +#+CALL: generate_f_interface(table=qmckl_slagel_splitting_args,rettyp=get_value("FRetType"),fname="qmckl_slagel_splitting_doc") + +#+RESULTS: +#+begin_src f90 :tangle (eval fh_func) :comments org :exports none +interface + integer(c_int32_t) function qmckl_slagel_splitting_doc & + (context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, & + Slater_inv, later_updates, later_index, later, determinant) & + bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + + integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: LDS + integer (c_int64_t) , intent(in) , value :: Dim + integer (c_int64_t) , intent(in) , value :: N_updates + real (c_double ) , intent(in) :: Updates(N_updates*LDS) + integer (c_int64_t) , intent(in) :: Updates_index(N_updates) + real (c_double ) , intent(in) , value :: breakdown + real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS) + real (c_double ) , intent(inout) :: later_updates(N_updates*LDS) + integer (c_int64_t) , intent(inout) :: later_index(N_updates) + integer (c_int64_t) , intent(inout) :: later + real (c_double ) , intent(inout) :: determinant + + end function qmckl_slagel_splitting_doc +end interface +#+end_src + +#+CALL: generate_f_interface(table=qmckl_slagel_splitting_args,rettyp=get_value("FRetType"),fname="qmckl_slagel_splitting") + +#+RESULTS: +#+begin_src f90 :tangle (eval fh_func) :comments org :exports none +interface + integer(c_int32_t) function qmckl_slagel_splitting & + (context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, & + Slater_inv, later_updates, later_index, later, determinant) & + bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + + integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: LDS + integer (c_int64_t) , intent(in) , value :: Dim + integer (c_int64_t) , intent(in) , value :: N_updates + real (c_double ) , intent(in) :: Updates(N_updates*LDS) + integer (c_int64_t) , intent(in) :: Updates_index(N_updates) + real (c_double ) , intent(in) , value :: breakdown + real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS) + real (c_double ) , intent(inout) :: later_updates(N_updates*LDS) + integer (c_int64_t) , intent(inout) :: later_index(N_updates) + integer (c_int64_t) , intent(inout) :: later + real (c_double ) , intent(inout) :: determinant + + end function qmckl_slagel_splitting +end interface +#+end_src + +*** Performance +This function cannot be used by itself and is used in Sherman-Morrison with update splitting and Woodbury 3x3 and 2x2 +with Sherman-Morrison and update splitting. Please look at the performance reccomendations for those two kernels. + + +* Sherman-Morrison with Slagel Splitting ** ~qmckl_sherman_morrison_splitting~ :PROPERTIES: :Name: qmckl_sherman_morrison_splitting @@ -914,37 +1511,76 @@ qmckl_exit_code qmckl_sherman_morrison_splitting_doc ( #+end_src *** C source -#+begin_src c :tangle (eval c) :comments org -qmckl_exit_code qmckl_sherman_morrison_splitting_hpc(const qmckl_context context, - const uint64_t LDS, - const uint64_t Dim, - const uint64_t N_updates, - const double* Updates, - const uint64_t* Updates_index, - const double breakdown, - double* Slater_inv, - double* determinant) { +#+NAME:splitting_switch-case_generator +#+begin_src python :noweb yes +text=""" +case {Dim}: { + rc = qmckl_slagel_splitting_{Dim}( + context, + N_updates, + Updates, + Updates_index, + breakdown, + Slater_inv, + later_updates, + later_index, &later, determinant); + break; +} +""" +result = [] +for Dim in <>: + Dim=str(Dim) + result.append(text.replace("{Dim}",Dim) ) + +return '\n'.join(result) +#+end_src + +#+begin_src c :tangle (eval c) :comments org :noweb yes +qmckl_exit_code qmckl_sherman_morrison_splitting_hpc( + const qmckl_context context, + const uint64_t LDS, + const uint64_t Dim, + const uint64_t N_updates, + const double* Updates, + const uint64_t* Updates_index, + const double breakdown, + double* Slater_inv, + double* determinant) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return qmckl_failwith(context, - QMCKL_NULL_CONTEXT, - "qmckl_sherman_morrison_splitting", - NULL); + return qmckl_failwith( + context, + QMCKL_NULL_CONTEXT, + "qmckl_sherman_morrison_splitting_hpc", + NULL); } double __attribute__((aligned(8))) later_updates[LDS * N_updates]; uint64_t later_index[N_updates]; uint64_t later = 0; - qmckl_exit_code rc = qmckl_slagel_splitting( - LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, - later_updates, later_index, &later, determinant); + qmckl_exit_code rc; + if (LDS == (1+(Dim-1)/SIMD_LENGTH)*SIMD_LENGTH) { + switch (Dim) { + <> + default: { + assert(0 == 1 && "TEMPLATE NOT IMPLEMENTED!"); + break; + } + } + } else { + rc = qmckl_slagel_splitting_hpc( + context, LDS, Dim, N_updates, Updates, Updates_index, + breakdown, Slater_inv, later_updates, + later_index, &later, determinant); + } if (rc != QMCKL_SUCCESS) return QMCKL_FAILURE; if (later > 0) { - qmckl_exit_code rc = qmckl_sherman_morrison_splitting( - context, LDS, Dim, later, later_updates, later_index, breakdown, - Slater_inv, determinant); + qmckl_exit_code rc = qmckl_sherman_morrison_splitting_hpc( + context, LDS, Dim, later, + later_updates, later_index, + breakdown, Slater_inv, determinant); if (rc != QMCKL_SUCCESS) return QMCKL_FAILURE; } @@ -952,56 +1588,60 @@ qmckl_exit_code qmckl_sherman_morrison_splitting_hpc(const qmckl_context context } #+end_src -#+begin_src c :tangle (eval c) :comments org -qmckl_exit_code qmckl_sherman_morrison_splitting(const qmckl_context context, - const uint64_t LDS, - const uint64_t Dim, - const uint64_t N_updates, - const double* Updates, - const uint64_t* Updates_index, - const double breakdown, - double* Slater_inv, - double* determinant) { +#+begin_src c :tangle (eval c) :comment org +qmckl_exit_code qmckl_sherman_morrison_splitting( + const qmckl_context context, + const uint64_t LDS, + const uint64_t Dim, + const uint64_t N_updates, + const double* Updates, + const uint64_t* Updates_index, + const double breakdown, + double* Slater_inv, + double* determinant) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return qmckl_failwith(context, - QMCKL_NULL_CONTEXT, - "qmckl_sherman_morrison_splitting", - NULL); + return qmckl_failwith( + context, + QMCKL_NULL_CONTEXT, + "qmckl_sherman_morrison_splitting", + NULL); } - #ifdef HAS_HPC - return qmckl_sherman_morrison_splitting_hpc(context, - LDS, - Dim, - N_updates, - Updates, - Updates_index, - breakdown, - Slater_inv, - determinant); + return qmckl_sherman_morrison_splitting_hpc( + context, + LDS, + Dim, + N_updates, + Updates, + Updates_index, + breakdown, + Slater_inv, + determinant); #else - // return qmckl_sherman_morrison_splitting_doc(context, - // LDS, - // Dim, - // N_updates, - // Updates, - // Updates_index, - // breakdown, - // Slater_inv, - // determinant); - return qmckl_sherman_morrison_splitting_hpc(context, - LDS, - Dim, - N_updates, - Updates, - Updates_index, - breakdown, - Slater_inv, - determinant); - #endif - - return QMCKL_SUCCESS; + // return qmckl_sherman_morrison_splitting_doc( + // context, + // LDS, + // Dim, + // N_updates, + // Updates, + // Updates_index, + // breakdown, + // Slater_inv, + // determinant); + return qmckl_sherman_morrison_splitting_hpc( + context, + LDS, + Dim, + N_updates, + Updates, + Updates_index, + breakdown, + Slater_inv, + determinant); + #endif + + return QMCKL_SUCCESS; } #+end_src @@ -1124,533 +1764,6 @@ assert(rc == QMCKL_SUCCESS); #+end_src -* Slagel Splitting -** ~qmckl_slagel_splitting~ -:PROPERTIES: -:Name: qmckl_slagel_splitting -:CRetType: qmckl_exit_code -:FRetType: qmckl_exit_code -:END: - -*** Introduction -~qmckl_slagel_splitting~ is the non-recursive, inner part of the 'Sherman-Morrison with update splitting'-kernel. -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}. - -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. - -*** API -#+NAME: qmckl_slagel_splitting_args -| Variable | Type | In/Out | Description | -|-----------------+-------------------------+--------+---------------------------------------------------------------| -| ~LDS~ | ~uint64_t~ | in | Leading dimension of Slater_inv | -| ~Dim~ | ~uint64_t~ | in | Dimension of Slater_inv | -| ~N_updates~ | ~uint64_t~ | in | Number of rank-1 updates to be applied to Slater_inv | -| ~Updates~ | ~double[N_updates*LDS]~ | in | Array containing the rank-1 updates | -| ~Updates_index~ | ~uint64_t[N_updates]~ | in | Array containing positions of the rank-1 updates | -| ~breakdown~ | ~double~ | in | Break-down parameter on which to fail or not | -| ~Slater_inv~ | ~double[Dim*LDS]~ | inout | Array containing the inverse Slater-matrix | -| ~later_updates~ | ~double[N_updates*LDS]~ | inout | Array containing the split updates for later | -| ~later_index~ | ~uint64_t[N_updates]~ | inout | Array containing the positions of the split updates for later | -| ~later~ | ~uint64_t~ | inout | Number of split updates for later | -| ~determinant~ | ~double~ | inout | Determinant of the Slater-matrix | - -*** Requirements - * ~LDS >= 2~ - * ~Dim >= 2~ - * ~N_updates >= 1~ - * ~Updates~ is allocated with $N_updates \times Dim$ elements - * ~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 - * ~later_updates~ is allocated with $later \times Dim$ elements - * ~later_index~ is allocated with $N_updates$ elements - * ~later >= 0~ - -*** Pedagogical kernel source (in Fortran) -The following source code written in Fortran is inteded to illustrate how the kernel works. Even though the kernel is -able to do numerically correct computations, it does not do it in the most efficient way possible. It should therefore -not be used in real workloads. - -#+begin_src f90 :tangle (eval f) -integer function qmckl_slagel_splitting_doc_f( & - lds, dim, & - nupdates, & - upds, & - updates_index, & - breakdown, & - s_inv, & - later_updates, & - later_index, & - later, & - determinant) result(info) - - use qmckl - implicit none - integer*8 , intent(in) :: lds, dim - integer*8 , intent(in) :: nupdates - integer*8 , intent(in) :: updates_index(nupdates) - real*8 , intent(in) :: upds(nupdates * lds) - real*8 , intent(in) :: breakdown - real*8 , intent(inout) :: s_inv(dim * lds) - real*8 , intent(inout) :: later_updates(nupdates * lds) - integer*8 , intent(inout) :: later_index(nupdates) - integer*8 , intent(inout) :: later - real*8 , intent(inout) :: determinant - - real*8 , dimension(lds, nupdates) :: Updates - real*8 , dimension(dim, lds) :: Inverse - - info = QMCKL_FAILURE - - ! Convert 'upds' and 's_inv' into the more easily readable Fortran - ! matrices 'Updates' and 'Inverse'. - call convert(upds, s_inv, Updates, Inverse, nupdates, lds, dim) - - ! YET TO BE IMPLEMENTED - - ! Copy updated inverse back to s_inv - call copy_back(Inverse, s_inv, lds, dim) - - info = QMCKL_SUCCESS - -end function qmckl_slagel_splitting_doc_f -#+end_src - -**** C interface to the pedagogical kernel (not directly exposed) -The following Fortran function ~qmckl_slagel_splitting_doc~ makes sure -that the pedagogical kernel ~qmckl_slagel_splitting_doc_f~, written in -Fortran, can be called from C using the ~ISO_C_BINDING~. The Fortran function -~qmckl_slagel_splitting_doc~ will be exposed in the header file 'qmckl.h' -for C users and in the module file 'qmckl_f.F90' for Fortran users. - -#+CALL: generate_c_interface(table=qmckl_slagel_splitting_args,rettyp=get_value("CRetType"),fname="qmckl_slagel_splitting_doc") - -#+RESULTS: -#+begin_src f90 :tangle (eval f) :comments org :exports none -integer(c_int32_t) function qmckl_slagel_splitting_doc & - (LDS, & - Dim, & - N_updates, & - Updates, & - Updates_index, & - breakdown, & - Slater_inv, & - later_updates, & - later_index, & - later, & - determinant) & - bind(C) result(info) - - use, intrinsic :: iso_c_binding - implicit none - - integer (c_int64_t) , intent(in) , value :: LDS - integer (c_int64_t) , intent(in) , value :: Dim - integer (c_int64_t) , intent(in) , value :: N_updates - real (c_double ) , intent(in) :: Updates(N_updates*LDS) - integer (c_int64_t) , intent(in) :: Updates_index(N_updates) - real (c_double ) , intent(in) , value :: breakdown - real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS) - real (c_double ) , intent(inout) :: later_updates(N_updates*LDS) - integer (c_int64_t) , intent(inout) :: later_index(N_updates) - integer (c_int64_t) , intent(inout) :: later - real (c_double ) , intent(inout) :: determinant - - integer(c_int32_t), external :: qmckl_slagel_splitting_doc_f - info = qmckl_slagel_splitting_doc_f & - (LDS, & - Dim, & - N_updates, & - Updates, & - Updates_index, & - breakdown, & - Slater_inv, & - later_updates, & - later_index, & - later, & - determinant) - -end function qmckl_slagel_splitting_doc -#+end_src - -*** C headers (exposed in qmckl.h) -#+CALL: generate_c_header(table=qmckl_slagel_splitting_args,rettyp=get_value("CRetType"),fname=get_value("Name")) - -#+RESULTS: -#+begin_src c :tangle (eval h_func) :comments org -qmckl_exit_code qmckl_slagel_splitting ( - const uint64_t LDS, - const uint64_t Dim, - const uint64_t N_updates, - const double* Updates, - const uint64_t* Updates_index, - const double breakdown, - double* Slater_inv, - double* later_updates, - uint64_t* later_index, - uint64_t* later, - double* determinant ); -#+end_src - -*** C sources -#+begin_src c :tangle (eval c) :comments org -qmckl_exit_code qmckl_slagel_splitting_hpc( - uint64_t LDS, - uint64_t Dim, - uint64_t N_updates, - const double* __restrict Updates, - const uint64_t* __restrict Updates_index, - const double breakdown, - double* __restrict Slater_inv, - double* __restrict later_updates, - uint64_t* __restrict later_index, - uint64_t* __restrict later, - double* __restrict determinant) { - - double __attribute__((aligned(8))) C[LDS]; - double __attribute__((aligned(8))) D[LDS]; - - uint64_t l = 0; - // For each update - while (l < N_updates) { - // C = S^{-1} x U_l - for (uint64_t i = 0; i < Dim; i++) { - C[i] = 0.0f; - IVDEP - ALIGNED - for (uint64_t j = 0; j < LDS; j++) { - C[i] += Slater_inv[i * LDS + j] * Updates[l * LDS + j]; - } - } - - // Denominator - const int cui = Updates_index[l] - 1; - double den = 1.0f + C[cui]; - if (fabs(den) < breakdown) { - // U_l = U_l / 2: split the update in 2 equal halves and save the - // second halve in later_updates - IVDEP - ALIGNED - for (uint64_t i = 0; i < LDS; i++) { - later_updates[*later * LDS + i] = Updates[l * LDS + i] * 0.5f; - C[i] *= 0.5f; - } - later_index[*later] = Updates_index[l]; - (*later)++; - - den = 1.0f + C[cui]; - } // From here onwards we continue with applying the first halve of the - // update to Slater_inv - double iden = 1.0f / den; - - if (determinant) - *determinant *= den; - - // D = v^T x S^{-1} : 1 x LDS - IVDEP - ALIGNED - for (uint64_t j = 0; j < LDS; j++) { - D[j] = Slater_inv[cui * LDS + j]; - } - - // S^{-1} = S^{-1} - C x D / den - for (uint64_t i = 0; i < Dim; i++) { - IVDEP - ALIGNED - for (uint64_t j = 0; j < LDS; j++) { - const double update = C[i] * D[j] * iden; - Slater_inv[i * LDS + j] -= update; - } - } - l += 1; - } - - return QMCKL_SUCCESS; -} - #+end_src - -#+NAME:slagel_splitting_template_code -#+begin_src c -static inline qmckl_exit_code qmckl_slagel_splitting_{Dim}( - uint64_t N_updates, - const double* __restrict Updates, - const uint64_t* __restrict Updates_index, - const double breakdown, - double* __restrict Slater_inv, - double* __restrict later_updates, - uint64_t* __restrict later_index, - uint64_t* __restrict later, - double* __restrict determinant) { - - double __attribute__((aligned(8))) C[D{Dim}_P]; - double __attribute__((aligned(8))) D[D{Dim}_P]; - - uint64_t l = 0; - // For each update - while (l < N_updates) { - // C = S^{-1} x U_l - for (uint64_t i = 0; i < {Dim}; i++) { - C[i] = 0.0f; - IVDEP - ALIGNED - for (uint64_t j = 0; j < D{Dim}_P; j++) { - C[i] += Slater_inv[i * D{Dim}_P + j] * Updates[l * D{Dim}_P + j]; - } - } - - // Denominator - const int cui = Updates_index[l] - 1; - double den = 1.0f + C[cui]; - if (fabs(den) < breakdown) { - // U_l = U_l / 2: split the update in 2 equal halves and save the - // second halve in later_updates - IVDEP - ALIGNED - for (uint64_t i = 0; i < D{Dim}_P; i++) { - later_updates[*later * D{Dim}_P + i] = Updates[l * D{Dim}_P + i] * 0.5f; - C[i] *= 0.5f; - } - later_index[*later] = Updates_index[l]; - (*later)++; - - den = 1.0f + C[cui]; - } // From here onwards we continue with applying the first halve of the - // update to Slater_inv - double iden = 1.0f / den; - - if (determinant) - *determinant *= den; - - // D = v^T x S^{-1} : 1 x D{Dim}_P - IVDEP - ALIGNED - for (uint64_t j = 0; j < D{Dim}_P; j++) { - D[j] = Slater_inv[cui * D{Dim}_P + j]; - } - - // S^{-1} = S^{-1} - C x D / den - for (uint64_t i = 0; i < {Dim}; i++) { - IVDEP - ALIGNED - for (uint64_t j = 0; j < D{Dim}_P; j++) { - const double update = C[i] * D[j] * iden; - Slater_inv[i * D{Dim}_P + j] -= update; - } - } - l += 1; - } - - return QMCKL_SUCCESS; -} -#+end_src - -#+NAME:slagel_splitting_kernel_generator -#+begin_src python :noweb yes -text=""" -<> -""" -result = [] -for Dim in <>: - Dim=str(Dim) - result.append(text.replace("{Dim}",Dim) ) - -return '\n'.join(result) -#+end_src - -#+NAME:slagel_splitting_switch-case_generator -#+begin_src python :noweb yes -text=""" -case {Dim}: - return qmckl_slagel_splitting_{Dim}( - N_updates, - Updates, - Updates_index, - breakdown, - Slater_inv, - later_updates, - later_index, - later, - determinant); -""" -result = [] -for Dim in <>: - Dim=str(Dim) - result.append(text.replace("{Dim}",Dim) ) - -return '\n'.join(result) -#+end_src - -#+begin_src c :tangle (eval c) :comments org :noweb yes -<> -#+end_src - -#+begin_src c :tangle (eval c) :comments org :noweb yes -qmckl_exit_code qmckl_slagel_splitting( - const uint64_t LDS, - const uint64_t Dim, - const uint64_t N_updates, - const double* Updates, - const uint64_t* Updates_index, - const double breakdown, - double* Slater_inv, - double* later_updates, - uint64_t* later_index, - uint64_t* later, - double* determinant) { - - #ifdef HAVE_HPC - if (LDS == (1+(Dim-1)/SIMD_LENGTH)*SIMD_LENGTH) { // Most cases - switch (Dim) { - <> - } - } - else { // Updating smaller sub-matrix - return qmckl_slagel_splitting_hpc( - LDS, - Dim, - N_updates, - Updates, - Updates_index, - breakdown, - Slater_inv, - later_updates, - later_index, - later, - determinant); - } - #else - // return qmckl_slagel_splitting_doc( - // LDS, - // Dim, - // N_updates, - // Updates, - // Updates_index, - // breakdown, - // Slater_inv, - // later_updates, - // later_index, - // later, - // determinant); - return qmckl_slagel_splitting_hpc( - LDS, - Dim, - N_updates, - Updates, - Updates_index, - breakdown, - Slater_inv, - later_updates, - later_index, - later, - determinant); - #endif - - return QMCKL_FAILURE; -} -#+end_src - -*** Fortran interfaces (exposed in qmckl_f.F90) -:PROPERTIES: -:Name: qmckl_slagel_splitting -:CRetType: qmckl_exit_code -:FRetType: qmckl_exit_code -:END: - -#+CALL: generate_f_interface(table=qmckl_slagel_splitting_args,rettyp=get_value("FRetType"),fname="qmckl_slagel_splitting_hpc") - -#+RESULTS: -#+begin_src f90 :tangle (eval fh_func) :comments org :exports none -interface - integer(c_int32_t) function qmckl_slagel_splitting_hpc & - (LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, later_updates, later_index, later, determinant) & - bind(C) - use, intrinsic :: iso_c_binding - import - implicit none - - integer (c_int64_t) , intent(in) , value :: LDS - integer (c_int64_t) , intent(in) , value :: Dim - integer (c_int64_t) , intent(in) , value :: N_updates - real (c_double ) , intent(in) :: Updates(N_updates*LDS) - integer (c_int64_t) , intent(in) :: Updates_index(N_updates) - real (c_double ) , intent(in) , value :: breakdown - real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS) - real (c_double ) , intent(inout) :: later_updates(N_updates*LDS) - integer (c_int64_t) , intent(inout) :: later_index(N_updates) - integer (c_int64_t) , intent(inout) :: later - real (c_double ) , intent(inout) :: determinant - - end function qmckl_slagel_splitting_hpc -end interface -#+end_src - -#+CALL: generate_f_interface(table=qmckl_slagel_splitting_args,rettyp=get_value("FRetType"),fname="qmckl_slagel_splitting_doc") - -#+RESULTS: -#+begin_src f90 :tangle (eval fh_func) :comments org :exports none -interface - integer(c_int32_t) function qmckl_slagel_splitting_doc & - (LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, later_updates, later_index, later, determinant) & - bind(C) - use, intrinsic :: iso_c_binding - import - implicit none - - integer (c_int64_t) , intent(in) , value :: LDS - integer (c_int64_t) , intent(in) , value :: Dim - integer (c_int64_t) , intent(in) , value :: N_updates - real (c_double ) , intent(in) :: Updates(N_updates*LDS) - integer (c_int64_t) , intent(in) :: Updates_index(N_updates) - real (c_double ) , intent(in) , value :: breakdown - real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS) - real (c_double ) , intent(inout) :: later_updates(N_updates*LDS) - integer (c_int64_t) , intent(inout) :: later_index(N_updates) - integer (c_int64_t) , intent(inout) :: later - real (c_double ) , intent(inout) :: determinant - - end function qmckl_slagel_splitting_doc -end interface -#+end_src - -#+CALL: generate_f_interface(table=qmckl_slagel_splitting_args,rettyp=get_value("FRetType"),fname="qmckl_slagel_splitting") - -#+RESULTS: -#+begin_src f90 :tangle (eval fh_func) :comments org :exports none -interface - integer(c_int32_t) function qmckl_slagel_splitting & - (LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, later_updates, later_index, later, determinant) & - bind(C) - use, intrinsic :: iso_c_binding - import - implicit none - - integer (c_int64_t) , intent(in) , value :: LDS - integer (c_int64_t) , intent(in) , value :: Dim - integer (c_int64_t) , intent(in) , value :: N_updates - real (c_double ) , intent(in) :: Updates(N_updates*LDS) - integer (c_int64_t) , intent(in) :: Updates_index(N_updates) - real (c_double ) , intent(in) , value :: breakdown - real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS) - real (c_double ) , intent(inout) :: later_updates(N_updates*LDS) - integer (c_int64_t) , intent(inout) :: later_index(N_updates) - integer (c_int64_t) , intent(inout) :: later - real (c_double ) , intent(inout) :: determinant - - end function qmckl_slagel_splitting -end interface -#+end_src - -*** Performance -This function cannot be used by itself and is used in Sherman-Morrison with update splitting and Woodbury 3x3 and 2x2 -with Sherman-Morrison and update splitting. Please look at the performance reccomendations for those two kernels. - * End of files From 8ba882675e3ee098b8cf949b8083b0519adf49da Mon Sep 17 00:00:00 2001 From: Francois Coppens Date: Mon, 20 Feb 2023 18:51:20 +0100 Subject: [PATCH 2/8] Renamed function prefixes. --- org/qmckl_sherman_morrison_woodbury.org | 244 ++++++++++++------------ 1 file changed, 122 insertions(+), 122 deletions(-) diff --git a/org/qmckl_sherman_morrison_woodbury.org b/org/qmckl_sherman_morrison_woodbury.org index 4c69d63..c91e185 100644 --- a/org/qmckl_sherman_morrison_woodbury.org +++ b/org/qmckl_sherman_morrison_woodbury.org @@ -34,9 +34,9 @@ range(2, 3) * Naïve Sherman-Morrison -** ~qmckl_sherman_morrison_naive~ +** ~qmckl_sm_naive~ :PROPERTIES: -:Name: qmckl_sherman_morrison_naive +:Name: qmckl_sm_naive :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: @@ -75,7 +75,7 @@ If the determinant of the Slater-matrix is passed, it will be updated to the det from applying the updates to the original matrix. *** API -#+NAME: qmckl_sherman_morrison_naive_args +#+NAME: qmckl_sm_naive_args | Variable | Type | In/Out | Description | |-----------------+-------------------------+--------+------------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | @@ -150,7 +150,7 @@ end subroutine copy_back #+end_src #+begin_src f90 :tangle (eval f) -integer function qmckl_sherman_morrison_naive_doc_f(context, & +integer function qmckl_sm_naive_doc_f(context, & lds, dim, & nupdates, & upds, & @@ -230,20 +230,20 @@ integer function qmckl_sherman_morrison_naive_doc_f(context, & info = QMCKL_SUCCESS -end function qmckl_sherman_morrison_naive_doc_f +end function qmckl_sm_naive_doc_f #+end_src **** C interface to the pedagogical kernel (not directly exposed) -The following Fortran function ~qmckl_sherman_morrison_naive_doc~ makes sure -that the pedagogical kernel ~qmckl_sherman_morrison_naive_doc_f~, written in -Fortran, can be called from C using the ~ISO_C_BINDING~. The Fortran function ~qmckl_sherman_morrison_naive_doc~ will be exposed in the header file 'qmckl.h' +The following Fortran function ~qmckl_sm_naive_doc~ makes sure +that the pedagogical kernel ~qmckl_sm_naive_doc_f~, written in +Fortran, can be called from C using the ~ISO_C_BINDING~. The Fortran function ~qmckl_sm_naive_doc~ will be exposed in the header file 'qmckl.h' for C users and in the module file 'qmckl_f.F90' for Fortran users. -#+CALL: generate_c_interface(table=qmckl_sherman_morrison_naive_args,rettyp=get_value("CRetType"),fname="qmckl_sherman_morrison_naive_doc") +#+CALL: generate_c_interface(table=qmckl_sm_naive_args,rettyp=get_value("CRetType"),fname="qmckl_sm_naive_doc") #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none -integer(c_int32_t) function qmckl_sherman_morrison_naive_doc & +integer(c_int32_t) function qmckl_sm_naive_doc & (context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) & bind(C) result(info) @@ -260,20 +260,20 @@ integer(c_int32_t) function qmckl_sherman_morrison_naive_doc & real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS) real (c_double ) , intent(inout) :: determinant - integer(c_int32_t), external :: qmckl_sherman_morrison_naive_doc_f - info = qmckl_sherman_morrison_naive_doc_f & + integer(c_int32_t), external :: qmckl_sm_naive_doc_f + info = qmckl_sm_naive_doc_f & (context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) -end function qmckl_sherman_morrison_naive_doc +end function qmckl_sm_naive_doc #+end_src *** C headers (exposed in qmckl.h) -#+CALL: generate_c_header(table=qmckl_sherman_morrison_naive_args,rettyp=get_value("CRetType"),fname=get_value("Name")) +#+CALL: generate_c_header(table=qmckl_sm_naive_args,rettyp=get_value("CRetType"),fname=get_value("Name")) #+RESULTS: #+begin_src c :tangle (eval h_func) : comments org -qmckl_exit_code qmckl_sherman_morrison_naive ( +qmckl_exit_code qmckl_sm_naive ( const qmckl_context context, const uint64_t LDS, const uint64_t Dim, @@ -285,11 +285,11 @@ qmckl_exit_code qmckl_sherman_morrison_naive ( double* determinant ); #+end_src -#+CALL: generate_c_header(table=qmckl_sherman_morrison_naive_args,rettyp=get_value("CRetType"),fname="qmckl_sherman_morrison_naive_hpc") +#+CALL: generate_c_header(table=qmckl_sm_naive_args,rettyp=get_value("CRetType"),fname="qmckl_sm_naive_hpc") #+RESULTS: #+begin_src c :tangle (eval h_private_func) :comments org -qmckl_exit_code qmckl_sherman_morrison_naive_hpc ( +qmckl_exit_code qmckl_sm_naive_hpc ( const qmckl_context context, const uint64_t LDS, const uint64_t Dim, @@ -301,11 +301,11 @@ qmckl_exit_code qmckl_sherman_morrison_naive_hpc ( double* determinant ); #+end_src -#+CALL: generate_c_header(table=qmckl_sherman_morrison_naive_args,rettyp=get_value("CRetType"),fname="qmckl_sherman_morrison_naive_doc") +#+CALL: generate_c_header(table=qmckl_sm_naive_args,rettyp=get_value("CRetType"),fname="qmckl_sm_naive_doc") #+RESULTS: #+begin_src c :tangle (eval h_func) :comments org -qmckl_exit_code qmckl_sherman_morrison_naive_doc ( +qmckl_exit_code qmckl_sm_naive_doc ( const qmckl_context context, const uint64_t LDS, const uint64_t Dim, @@ -344,8 +344,8 @@ Common includes and macros used by all the Sherman-Morrison-Woodbury kernels. #endif #+end_src -~qmckl_sherman_morrison_naive_hpc~ is a high performance variation of -~qmckl_sherman_morrison_naive~ written in C. It is used in cases when ~Dim~ is +~qmckl_sm_naive_hpc~ is a high performance variation of +~qmckl_sm_naive~ written in C. It is used in cases when ~Dim~ is smaller than the leading dimension ~LDS~, irrespective of whetether ~LDS~ includes zero padding to benefit from SIMD instructions or not. Cases like this include situations where one wants to apply updates to a square submatrix of the @@ -354,7 +354,7 @@ It takes advantage of memory aligned data and assumes no data dependencies inside the loops. The loops are fully vectorised whenever ~Dim~ is an integer multiple of ~SIMD_LEGTH~. #+begin_src c :tangle (eval c) :comments org -qmckl_exit_code qmckl_sherman_morrison_naive_hpc( +qmckl_exit_code qmckl_sm_naive_hpc( const qmckl_context context, const uint64_t LDS, const uint64_t Dim, @@ -368,7 +368,7 @@ qmckl_exit_code qmckl_sherman_morrison_naive_hpc( if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith( context, QMCKL_NULL_CONTEXT, - "qmckl_sherman_morrison_naive_hpc", + "qmckl_sm_naive_hpc", NULL); } @@ -423,10 +423,10 @@ qmckl_exit_code qmckl_sherman_morrison_naive_hpc( } #+end_src -~qmckl_exit_code qmckl_sherman_morrison_naive_{Dim}~ is a C function-template that is used to genereate instances of C fucntions based on the range given above. The advantage of this method is that for each of these instances all the dimensions and loop-bounds are known at compile time, allowing the compiler to optimize more aggressively. +~qmckl_exit_code qmckl_sm_naive_{Dim}~ is a C function-template that is used to genereate instances of C fucntions based on the range given above. The advantage of this method is that for each of these instances all the dimensions and loop-bounds are known at compile time, allowing the compiler to optimize more aggressively. #+NAME:naive_template_code #+begin_src c - static inline qmckl_exit_code qmckl_sherman_morrison_naive_{Dim}( + static inline qmckl_exit_code qmckl_sm_naive_{Dim}( const qmckl_context context, const uint64_t N_updates, const double* __restrict Updates, @@ -438,7 +438,7 @@ qmckl_exit_code qmckl_sherman_morrison_naive_hpc( if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith(context, QMCKL_NULL_CONTEXT, - "qmckl_sherman_morrison_naive_{Dim}", + "qmckl_sm_naive_{Dim}", NULL); } @@ -516,7 +516,7 @@ Python script that generated C switch cases that call individual kernel instance #+begin_src python :noweb yes text=""" case {Dim}: - return qmckl_sherman_morrison_naive_{Dim}(context, + return qmckl_sm_naive_{Dim}(context, N_updates, Updates, Updates_index, @@ -536,9 +536,9 @@ return '\n'.join(result) <> #+end_src -~qmckl_sherman_morrison_naive~ is a generic function that contains decision making logic that calls the proper kernel based on the used library configuration (~--enable-doc~ and ~--enable-hpc~) and the passed array dimensions ~LDS~ and ~Dim~. +~qmckl_sm_naive~ is a generic function that contains decision making logic that calls the proper kernel based on the used library configuration (~--enable-doc~ and ~--enable-hpc~) and the passed array dimensions ~LDS~ and ~Dim~. #+begin_src c :tangle (eval c) :comments org :noweb yes -qmckl_exit_code qmckl_sherman_morrison_naive(const qmckl_context context, +qmckl_exit_code qmckl_sm_naive(const qmckl_context context, const uint64_t LDS, const uint64_t Dim, const uint64_t N_updates, @@ -551,7 +551,7 @@ qmckl_exit_code qmckl_sherman_morrison_naive(const qmckl_context context, if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith(context, QMCKL_NULL_CONTEXT, - "qmckl_sherman_morrison_naive", + "qmckl_sm_naive", NULL); } @@ -562,7 +562,7 @@ qmckl_exit_code qmckl_sherman_morrison_naive(const qmckl_context context, } } else { // Updating smaller sub-matrix - return qmckl_sherman_morrison_naive_hpc(context, + return qmckl_sm_naive_hpc(context, LDS, Dim, N_updates, @@ -573,7 +573,7 @@ qmckl_exit_code qmckl_sherman_morrison_naive(const qmckl_context context, determinant); } #else - return qmckl_sherman_morrison_naive_doc(context, + return qmckl_sm_naive_doc(context, LDS, Dim, N_updates, @@ -590,17 +590,17 @@ qmckl_exit_code qmckl_sherman_morrison_naive(const qmckl_context context, *** Fortran interfaces (exposed in qmckl_f.F90) :PROPERTIES: - :Name: qmckl_sherman_morrison_naive + :Name: qmckl_sm_naive :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: -#+CALL: generate_f_interface(table=qmckl_sherman_morrison_naive_args,rettyp=get_value("FRetType"),fname="qmckl_sherman_morrison_naive_hpc") +#+CALL: generate_f_interface(table=qmckl_sm_naive_args,rettyp=get_value("FRetType"),fname="qmckl_sm_naive_hpc") #+RESULTS: #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface - integer(c_int32_t) function qmckl_sherman_morrison_naive_hpc & + integer(c_int32_t) function qmckl_sm_naive_hpc & (context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) & bind(C) use, intrinsic :: iso_c_binding @@ -617,16 +617,16 @@ interface real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS) real (c_double ) , intent(inout) :: determinant - end function qmckl_sherman_morrison_naive_hpc + end function qmckl_sm_naive_hpc end interface #+end_src -#+CALL: generate_f_interface(table=qmckl_sherman_morrison_naive_args,rettyp=get_value("FRetType"),fname="qmckl_sherman_morrison_naive_doc") +#+CALL: generate_f_interface(table=qmckl_sm_naive_args,rettyp=get_value("FRetType"),fname="qmckl_sm_naive_doc") #+RESULTS: #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface - integer(c_int32_t) function qmckl_sherman_morrison_naive_doc & + integer(c_int32_t) function qmckl_sm_naive_doc & (context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) & bind(C) use, intrinsic :: iso_c_binding @@ -643,16 +643,16 @@ interface real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS) real (c_double ) , intent(inout) :: determinant - end function qmckl_sherman_morrison_naive_doc + end function qmckl_sm_naive_doc end interface #+end_src -#+CALL: generate_f_interface(table=qmckl_sherman_morrison_naive_args,rettyp=get_value("FRetType"),fname="qmckl_sherman_morrison_naive") +#+CALL: generate_f_interface(table=qmckl_sm_naive_args,rettyp=get_value("FRetType"),fname="qmckl_sm_naive") #+RESULTS: #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface - integer(c_int32_t) function qmckl_sherman_morrison_naive & + integer(c_int32_t) function qmckl_sm_naive & (context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) & bind(C) use, intrinsic :: iso_c_binding @@ -669,7 +669,7 @@ interface real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS) real (c_double ) , intent(inout) :: determinant - end function qmckl_sherman_morrison_naive + end function qmckl_sm_naive end interface #+end_src @@ -701,7 +701,7 @@ assert(Slater_inv1 != NULL); // original determinant of Slater1 (before applying updates) double det = 3.407025646103221e-10; -rc = qmckl_sherman_morrison_naive(context, +rc = qmckl_sm_naive(context, LDS, Dim, N_updates1, @@ -738,15 +738,15 @@ assert(rc == QMCKL_SUCCESS); * Sherman-Morrison with Slagel Splitting (core) -** ~qmckl_slagel_splitting~ +** ~qmckl_sm_splitting_core~ :PROPERTIES: -:Name: qmckl_slagel_splitting +:Name: qmckl_sm_splitting_core :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: *** Introduction -~qmckl_slagel_splitting~ is the inner core part of 'Sherman-Morrison with update splitting' in the next section. +~qmckl_sm_splitting_core~ is the inner core part of 'Sherman-Morrison with update splitting' in the next section. It is not normally used by itself but it is possible to use it nonetheless. It has three extra parameters in its API: @@ -754,14 +754,14 @@ It has three extra parameters in its API: - ~later_index~ initially empty array that will contain the row/column numbers of the updates that were split during execution - ~later~ initially zero integer that records the number of updates that were split during exection. -It is up to the user to decide what to do with these updates once the kernel returns. Normally ~qmckl_slagel_splitting~ is -used as the core part of a recursive function, as is done in ~qmckl_sherman_morrison_splitting~ or as part of a more complex +It is up to the user to decide what to do with these updates once the kernel returns. Normally ~qmckl_sm_splitting_core~ is +used as the core part of a recursive function, as is done in ~qmckl_sm_splitting~ or as part of a more complex kernel like ~qmckl_sherman_morrison_smw32s~. If the determinant is passed it will only be partially updated if there were any update splits. *** API -#+NAME: qmckl_slagel_splitting_args +#+NAME: qmckl_sm_splitting_core_args | Variable | Type | In/Out | Description | |-----------------+-------------------------+--------+---------------------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | @@ -795,7 +795,7 @@ able to do numerically correct computations, it does not do it in the most effic not be used in real workloads. #+begin_src f90 :tangle (eval f) -integer function qmckl_slagel_splitting_doc_f( & +integer function qmckl_sm_splitting_core_doc_f( & context, & lds, dim, & nupdates, & @@ -843,21 +843,21 @@ integer function qmckl_slagel_splitting_doc_f( & info = QMCKL_SUCCESS -end function qmckl_slagel_splitting_doc_f +end function qmckl_sm_splitting_core_doc_f #+end_src **** C interface to the pedagogical kernel (not directly exposed) -The following Fortran function ~qmckl_slagel_splitting_doc~ makes sure -that the pedagogical kernel ~qmckl_slagel_splitting_doc_f~, written in +The following Fortran function ~qmckl_sm_splitting_core_doc~ makes sure +that the pedagogical kernel ~qmckl_sm_splitting_core_doc_f~, written in Fortran, can be called from C using the ~ISO_C_BINDING~. The Fortran function -~qmckl_slagel_splitting_doc~ will be exposed in the header file 'qmckl.h' +~qmckl_sm_splitting_core_doc~ will be exposed in the header file 'qmckl.h' for C users and in the module file 'qmckl_f.F90' for Fortran users. -#+CALL: generate_c_interface(table=qmckl_slagel_splitting_args,rettyp=get_value("CRetType"),fname="qmckl_slagel_splitting_doc") +#+CALL: generate_c_interface(table=qmckl_sm_splitting_core_args,rettyp=get_value("CRetType"),fname="qmckl_sm_splitting_core_doc") #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none -integer(c_int32_t) function qmckl_slagel_splitting_doc & +integer(c_int32_t) function qmckl_sm_splitting_core_doc & (context, & LDS, & Dim, & @@ -888,8 +888,8 @@ integer(c_int32_t) function qmckl_slagel_splitting_doc & integer (c_int64_t) , intent(inout) :: later real (c_double ) , intent(inout) :: determinant - integer(c_int32_t), external :: qmckl_slagel_splitting_doc_f - info = qmckl_slagel_splitting_doc_f & + integer(c_int32_t), external :: qmckl_sm_splitting_core_doc_f + info = qmckl_sm_splitting_core_doc_f & (context, & LDS, & Dim, & @@ -903,15 +903,15 @@ integer(c_int32_t) function qmckl_slagel_splitting_doc & later, & determinant) -end function qmckl_slagel_splitting_doc +end function qmckl_sm_splitting_core_doc #+end_src *** C headers (exposed in qmckl.h) -#+CALL: generate_c_header(table=qmckl_slagel_splitting_args,rettyp=get_value("CRetType"),fname=get_value("Name")) +#+CALL: generate_c_header(table=qmckl_sm_splitting_core_args,rettyp=get_value("CRetType"),fname=get_value("Name")) #+RESULTS: #+begin_src c :tangle (eval h_func) :comments org -qmckl_exit_code qmckl_slagel_splitting ( +qmckl_exit_code qmckl_sm_splitting_core ( const qmckl_context context, const uint64_t LDS, const uint64_t Dim, @@ -926,11 +926,11 @@ qmckl_exit_code qmckl_slagel_splitting ( double* determinant ); #+end_src -#+CALL: generate_c_header(table=qmckl_slagel_splitting_args,rettyp=get_value("CRetType"),fname="qmckl_slagel_splitting_doc") +#+CALL: generate_c_header(table=qmckl_sm_splitting_core_args,rettyp=get_value("CRetType"),fname="qmckl_sm_splitting_core_doc") #+RESULTS: #+begin_src c :tangle (eval h_func) :comments org -qmckl_exit_code qmckl_slagel_splitting_doc ( +qmckl_exit_code qmckl_sm_splitting_core_doc ( const qmckl_context context, const uint64_t LDS, const uint64_t Dim, @@ -947,7 +947,7 @@ qmckl_exit_code qmckl_slagel_splitting_doc ( *** C sources #+begin_src c :tangle (eval c) :comments org -qmckl_exit_code qmckl_slagel_splitting_hpc( +qmckl_exit_code qmckl_sm_splitting_core_hpc( const qmckl_context context, uint64_t LDS, uint64_t Dim, @@ -965,7 +965,7 @@ qmckl_exit_code qmckl_slagel_splitting_hpc( return qmckl_failwith( context, QMCKL_NULL_CONTEXT, - "qmckl_slagel_splitting_hpc", + "qmckl_sm_splitting_core_hpc", NULL); } @@ -1033,7 +1033,7 @@ qmckl_exit_code qmckl_slagel_splitting_hpc( #+NAME:slagel_splitting_template_code #+begin_src c -static inline qmckl_exit_code qmckl_slagel_splitting_{Dim}( +static inline qmckl_exit_code qmckl_sm_splitting_core_{Dim}( const qmckl_context context, uint64_t N_updates, const double* __restrict Updates, @@ -1049,7 +1049,7 @@ static inline qmckl_exit_code qmckl_slagel_splitting_{Dim}( return qmckl_failwith( context, QMCKL_NULL_CONTEXT, - "qmckl_slagel_splitting_{Dim}", + "qmckl_sm_splitting_core_{Dim}", NULL); } @@ -1132,7 +1132,7 @@ return '\n'.join(result) #+begin_src python :noweb yes text=""" case {Dim}: { - return qmckl_slagel_splitting_{Dim}( + return qmckl_sm_splitting_core_{Dim}( context, N_updates, Updates, @@ -1159,7 +1159,7 @@ return '\n'.join(result) #+end_src #+begin_src c :tangle (eval c) :comments org :noweb yes -qmckl_exit_code qmckl_slagel_splitting( +qmckl_exit_code qmckl_sm_splitting_core( const qmckl_context context, const uint64_t LDS, const uint64_t Dim, @@ -1184,7 +1184,7 @@ qmckl_exit_code qmckl_slagel_splitting( } } else { // Updating smaller sub-matrix - return qmckl_slagel_splitting_hpc( + return qmckl_sm_splitting_core_hpc( context, LDS, Dim, @@ -1199,7 +1199,7 @@ qmckl_exit_code qmckl_slagel_splitting( determinant); } #else - // return qmckl_slagel_splitting_doc( + // return qmckl_sm_splitting_core_doc( // context, // LDS, // Dim, @@ -1212,7 +1212,7 @@ qmckl_exit_code qmckl_slagel_splitting( // later_index, // later, // determinant); - return qmckl_slagel_splitting_hpc( + return qmckl_sm_splitting_core_hpc( context, LDS, Dim, @@ -1233,17 +1233,17 @@ qmckl_exit_code qmckl_slagel_splitting( *** Fortran interfaces (exposed in qmckl_f.F90) :PROPERTIES: -:Name: qmckl_slagel_splitting +:Name: qmckl_sm_splitting_core :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: -#+CALL: generate_f_interface(table=qmckl_slagel_splitting_args,rettyp=get_value("FRetType"),fname="qmckl_slagel_splitting_hpc") +#+CALL: generate_f_interface(table=qmckl_sm_splitting_core_args,rettyp=get_value("FRetType"),fname="qmckl_sm_splitting_core_hpc") #+RESULTS: #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface - integer(c_int32_t) function qmckl_slagel_splitting_hpc & + integer(c_int32_t) function qmckl_sm_splitting_core_hpc & (context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, & Slater_inv, later_updates, later_index, later, determinant) & bind(C) @@ -1264,16 +1264,16 @@ interface integer (c_int64_t) , intent(inout) :: later real (c_double ) , intent(inout) :: determinant - end function qmckl_slagel_splitting_hpc + end function qmckl_sm_splitting_core_hpc end interface #+end_src -#+CALL: generate_f_interface(table=qmckl_slagel_splitting_args,rettyp=get_value("FRetType"),fname="qmckl_slagel_splitting_doc") +#+CALL: generate_f_interface(table=qmckl_sm_splitting_core_args,rettyp=get_value("FRetType"),fname="qmckl_sm_splitting_core_doc") #+RESULTS: #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface - integer(c_int32_t) function qmckl_slagel_splitting_doc & + integer(c_int32_t) function qmckl_sm_splitting_core_doc & (context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, & Slater_inv, later_updates, later_index, later, determinant) & bind(C) @@ -1294,16 +1294,16 @@ interface integer (c_int64_t) , intent(inout) :: later real (c_double ) , intent(inout) :: determinant - end function qmckl_slagel_splitting_doc + end function qmckl_sm_splitting_core_doc end interface #+end_src -#+CALL: generate_f_interface(table=qmckl_slagel_splitting_args,rettyp=get_value("FRetType"),fname="qmckl_slagel_splitting") +#+CALL: generate_f_interface(table=qmckl_sm_splitting_core_args,rettyp=get_value("FRetType"),fname="qmckl_sm_splitting_core") #+RESULTS: #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface - integer(c_int32_t) function qmckl_slagel_splitting & + integer(c_int32_t) function qmckl_sm_splitting_core & (context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, & Slater_inv, later_updates, later_index, later, determinant) & bind(C) @@ -1324,7 +1324,7 @@ interface integer (c_int64_t) , intent(inout) :: later real (c_double ) , intent(inout) :: determinant - end function qmckl_slagel_splitting + end function qmckl_sm_splitting_core end interface #+end_src @@ -1334,9 +1334,9 @@ with Sherman-Morrison and update splitting. Please look at the performance recco * Sherman-Morrison with Slagel Splitting -** ~qmckl_sherman_morrison_splitting~ +** ~qmckl_sm_splitting~ :PROPERTIES: -:Name: qmckl_sherman_morrison_splitting +:Name: qmckl_sm_splitting :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: @@ -1357,7 +1357,7 @@ If the determinant of the Slater-matrix is passed, it will be updated to the det from applying the updates to the original matrix. *** API -#+NAME: qmckl_sherman_morrison_splitting_args +#+NAME: qmckl_sm_splitting_args | Variable | Type | In/Out | Description | |---------------+-----------------------+--------+------------------------------------------------------| | context | qmckl_context | in | Global state | @@ -1386,7 +1386,7 @@ able to do numerically correct computations, it does not do it in the most effic not be used in real workloads. #+begin_src f90 :tangle (eval f) -integer function qmckl_sherman_morrison_splitting_doc_f(context, & +integer function qmckl_sm_splitting_doc_f(context, & lds, dim, & nupdates, & upds, & @@ -1422,21 +1422,21 @@ integer function qmckl_sherman_morrison_splitting_doc_f(context, & info = QMCKL_SUCCESS -end function qmckl_sherman_morrison_splitting_doc_f +end function qmckl_sm_splitting_doc_f #+end_src **** C interface to the pedagogical kernel (not directly exposed) -The following Fortran function ~qmckl_slagel_splitting_doc~ makes sure -that the pedagogical kernel ~qmckl_slagel_splitting_doc_f~, written in +The following Fortran function ~qmckl_sm_splitting_core_doc~ makes sure +that the pedagogical kernel ~qmckl_sm_splitting_core_doc_f~, written in Fortran, can be called from C using the ~ISO_C_BINDING~. The Fortran function -~qmckl_slagel_splitting_doc~ will be exposed in the header file 'qmckl.h' +~qmckl_sm_splitting_core_doc~ will be exposed in the header file 'qmckl.h' for C users and in the module file 'qmckl_f.F90' for Fortran users. -#+CALL: generate_c_interface(table=qmckl_sherman_morrison_splitting_args,rettyp=get_value("CRetType"),fname="qmckl_sherman_morrison_splitting_doc") +#+CALL: generate_c_interface(table=qmckl_sm_splitting_args,rettyp=get_value("CRetType"),fname="qmckl_sm_splitting_doc") #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none -integer(c_int32_t) function qmckl_sherman_morrison_splitting_doc & +integer(c_int32_t) function qmckl_sm_splitting_doc & (context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) & bind(C) result(info) @@ -1453,20 +1453,20 @@ integer(c_int32_t) function qmckl_sherman_morrison_splitting_doc & real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS) real (c_double ) , intent(inout) :: determinant - integer(c_int32_t), external :: qmckl_sherman_morrison_splitting_doc_f - info = qmckl_sherman_morrison_splitting_doc_f & + integer(c_int32_t), external :: qmckl_sm_splitting_doc_f + info = qmckl_sm_splitting_doc_f & (context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) -end function qmckl_sherman_morrison_splitting_doc +end function qmckl_sm_splitting_doc #+end_src *** C headers (exposed in qmckl.h) -#+CALL: generate_c_header(table=qmckl_sherman_morrison_splitting_args,rettyp=get_value("CRetType"),fname=get_value("Name")) +#+CALL: generate_c_header(table=qmckl_sm_splitting_args,rettyp=get_value("CRetType"),fname=get_value("Name")) #+RESULTS: #+begin_src c :tangle (eval h_func) :comments org -qmckl_exit_code qmckl_sherman_morrison_splitting ( +qmckl_exit_code qmckl_sm_splitting ( const qmckl_context context, const uint64_t LDS, const uint64_t Dim, @@ -1478,11 +1478,11 @@ qmckl_exit_code qmckl_sherman_morrison_splitting ( double* determinant ); #+end_src -#+CALL: generate_c_header(table=qmckl_sherman_morrison_splitting_args,rettyp=get_value("CRetType"),fname="qmckl_sherman_morrison_splitting_hpc") +#+CALL: generate_c_header(table=qmckl_sm_splitting_args,rettyp=get_value("CRetType"),fname="qmckl_sm_splitting_hpc") #+RESULTS: #+begin_src c :tangle (eval h_func) :comments org -qmckl_exit_code qmckl_sherman_morrison_splitting_hpc ( +qmckl_exit_code qmckl_sm_splitting_hpc ( const qmckl_context context, const uint64_t LDS, const uint64_t Dim, @@ -1494,11 +1494,11 @@ qmckl_exit_code qmckl_sherman_morrison_splitting_hpc ( double* determinant ); #+end_src -#+CALL: generate_c_header(table=qmckl_sherman_morrison_splitting_args,rettyp=get_value("CRetType"),fname="qmckl_sherman_morrison_splitting_doc") +#+CALL: generate_c_header(table=qmckl_sm_splitting_args,rettyp=get_value("CRetType"),fname="qmckl_sm_splitting_doc") #+RESULTS: #+begin_src c :tangle (eval h_func) :comments org -qmckl_exit_code qmckl_sherman_morrison_splitting_doc ( +qmckl_exit_code qmckl_sm_splitting_doc ( const qmckl_context context, const uint64_t LDS, const uint64_t Dim, @@ -1515,7 +1515,7 @@ qmckl_exit_code qmckl_sherman_morrison_splitting_doc ( #+begin_src python :noweb yes text=""" case {Dim}: { - rc = qmckl_slagel_splitting_{Dim}( + rc = qmckl_sm_splitting_core_{Dim}( context, N_updates, Updates, @@ -1536,7 +1536,7 @@ return '\n'.join(result) #+end_src #+begin_src c :tangle (eval c) :comments org :noweb yes -qmckl_exit_code qmckl_sherman_morrison_splitting_hpc( +qmckl_exit_code qmckl_sm_splitting_hpc( const qmckl_context context, const uint64_t LDS, const uint64_t Dim, @@ -1551,7 +1551,7 @@ qmckl_exit_code qmckl_sherman_morrison_splitting_hpc( return qmckl_failwith( context, QMCKL_NULL_CONTEXT, - "qmckl_sherman_morrison_splitting_hpc", + "qmckl_sm_splitting_hpc", NULL); } @@ -1569,7 +1569,7 @@ qmckl_exit_code qmckl_sherman_morrison_splitting_hpc( } } } else { - rc = qmckl_slagel_splitting_hpc( + rc = qmckl_sm_splitting_core_hpc( context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, later_updates, later_index, &later, determinant); @@ -1577,7 +1577,7 @@ qmckl_exit_code qmckl_sherman_morrison_splitting_hpc( if (rc != QMCKL_SUCCESS) return QMCKL_FAILURE; if (later > 0) { - qmckl_exit_code rc = qmckl_sherman_morrison_splitting_hpc( + qmckl_exit_code rc = qmckl_sm_splitting_hpc( context, LDS, Dim, later, later_updates, later_index, breakdown, Slater_inv, determinant); @@ -1589,7 +1589,7 @@ qmckl_exit_code qmckl_sherman_morrison_splitting_hpc( #+end_src #+begin_src c :tangle (eval c) :comment org -qmckl_exit_code qmckl_sherman_morrison_splitting( +qmckl_exit_code qmckl_sm_splitting( const qmckl_context context, const uint64_t LDS, const uint64_t Dim, @@ -1604,11 +1604,11 @@ qmckl_exit_code qmckl_sherman_morrison_splitting( return qmckl_failwith( context, QMCKL_NULL_CONTEXT, - "qmckl_sherman_morrison_splitting", + "qmckl_sm_splitting", NULL); } #ifdef HAS_HPC - return qmckl_sherman_morrison_splitting_hpc( + return qmckl_sm_splitting_hpc( context, LDS, Dim, @@ -1619,7 +1619,7 @@ qmckl_exit_code qmckl_sherman_morrison_splitting( Slater_inv, determinant); #else - // return qmckl_sherman_morrison_splitting_doc( + // return qmckl_sm_splitting_doc( // context, // LDS, // Dim, @@ -1629,7 +1629,7 @@ qmckl_exit_code qmckl_sherman_morrison_splitting( // breakdown, // Slater_inv, // determinant); - return qmckl_sherman_morrison_splitting_hpc( + return qmckl_sm_splitting_hpc( context, LDS, Dim, @@ -1647,17 +1647,17 @@ qmckl_exit_code qmckl_sherman_morrison_splitting( *** Fortran interfaces (exposed in qmckl_f.F90) :PROPERTIES: - :Name: qmckl_sherman_morrison_naive + :Name: qmckl_sm_naive :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: -#+CALL: generate_f_interface(table=qmckl_sherman_morrison_splitting_args,rettyp=get_value("FRetType"),fname="qmckl_sherman_morrison_splitting_hpc") +#+CALL: generate_f_interface(table=qmckl_sm_splitting_args,rettyp=get_value("FRetType"),fname="qmckl_sm_splitting_hpc") #+RESULTS: #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface - integer(c_int32_t) function qmckl_sherman_morrison_splitting_hpc & + integer(c_int32_t) function qmckl_sm_splitting_hpc & (context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) & bind(C) use, intrinsic :: iso_c_binding @@ -1674,16 +1674,16 @@ interface real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS) real (c_double ) , intent(inout) :: determinant - end function qmckl_sherman_morrison_splitting_hpc + end function qmckl_sm_splitting_hpc end interface #+end_src -#+CALL: generate_f_interface(table=qmckl_sherman_morrison_splitting_args,rettyp=get_value("FRetType"),fname="qmckl_sherman_morrison_splitting_doc") +#+CALL: generate_f_interface(table=qmckl_sm_splitting_args,rettyp=get_value("FRetType"),fname="qmckl_sm_splitting_doc") #+RESULTS: #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface - integer(c_int32_t) function qmckl_sherman_morrison_splitting_doc & + integer(c_int32_t) function qmckl_sm_splitting_doc & (context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) & bind(C) use, intrinsic :: iso_c_binding @@ -1700,16 +1700,16 @@ interface real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS) real (c_double ) , intent(inout) :: determinant - end function qmckl_sherman_morrison_splitting_doc + end function qmckl_sm_splitting_doc end interface #+end_src -#+CALL: generate_f_interface(table=qmckl_sherman_morrison_splitting_args,rettyp=get_value("FRetType"),fname="qmckl_sherman_morrison_splitting") +#+CALL: generate_f_interface(table=qmckl_sm_splitting_args,rettyp=get_value("FRetType"),fname="qmckl_sm_splitting") #+RESULTS: #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface - integer(c_int32_t) function qmckl_sherman_morrison_splitting & + integer(c_int32_t) function qmckl_sm_splitting & (context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) & bind(C) use, intrinsic :: iso_c_binding @@ -1726,7 +1726,7 @@ interface real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS) real (c_double ) , intent(inout) :: determinant - end function qmckl_sherman_morrison_splitting + end function qmckl_sm_splitting end interface #+end_src @@ -1739,7 +1739,7 @@ assert(Updates3 != NULL); assert(Updates_index3 != NULL); assert(Slater_inv3_2 != NULL); det = -1.23743195512859e-09; -rc = qmckl_sherman_morrison_splitting(context, LDS, Dim, N_updates3, Updates3, Updates_index3, breakdown, Slater_inv3_2, &det); +rc = qmckl_sm_splitting(context, LDS, 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 j = 0; j < Dim; j++) { From 8216f682b38ce0e7598fe4d2a0757656603de873 Mon Sep 17 00:00:00 2001 From: Francois Coppens Date: Tue, 21 Feb 2023 19:27:23 +0100 Subject: [PATCH 3/8] Strange Fortran type error... --- org/qmckl_sherman_morrison_woodbury.org | 363 ++++++++++++++---------- 1 file changed, 217 insertions(+), 146 deletions(-) diff --git a/org/qmckl_sherman_morrison_woodbury.org b/org/qmckl_sherman_morrison_woodbury.org index c91e185..d323cab 100644 --- a/org/qmckl_sherman_morrison_woodbury.org +++ b/org/qmckl_sherman_morrison_woodbury.org @@ -132,7 +132,7 @@ end subroutine convert #+end_src #+begin_src f90 :tangle (eval f) :comment org :exports none -subroutine copy_back(Inverse, s_inv, lds, dim) +subroutine copy_back_inv(Inverse, s_inv, lds, dim) implicit none integer*8 , intent(in) :: lds, dim real*8 , intent(in) , dimension(dim, lds) :: Inverse @@ -146,7 +146,25 @@ subroutine copy_back(Inverse, s_inv, lds, dim) s_inv((i - 1) * lds + j) = Inverse(i, j) end do end do -end subroutine copy_back +end subroutine copy_back_inv +#+end_src + +#+begin_src f90 :tangle (eval f) :comment org :exports none +subroutine copy_back_lu(Later_updates, later_upds, lds, nupdates) + implicit none + integer*8 , intent(in) :: lds, nupdates + real*8 , intent(in) , dimension(nupdates, lds) :: Later_updates + real*8 , intent(out) :: later_upds(nupdates * lds) + + integer*8 :: i, j + + ! Copy updated inverse back to s_inv + do i = 1, nupdates + do j = 1, lds + later_upds((i - 1) * lds + j) = Later_updates(i, j) + end do + end do +end subroutine copy_back_lu #+end_src #+begin_src f90 :tangle (eval f) @@ -539,20 +557,21 @@ return '\n'.join(result) ~qmckl_sm_naive~ is a generic function that contains decision making logic that calls the proper kernel based on the used library configuration (~--enable-doc~ and ~--enable-hpc~) and the passed array dimensions ~LDS~ and ~Dim~. #+begin_src c :tangle (eval c) :comments org :noweb yes qmckl_exit_code qmckl_sm_naive(const qmckl_context context, - const uint64_t LDS, - const uint64_t Dim, - const uint64_t N_updates, - const double* Updates, - const uint64_t* Updates_index, - const double breakdown, - double* Slater_inv, - double* determinant) { + const uint64_t LDS, + const uint64_t Dim, + const uint64_t N_updates, + const double* Updates, + const uint64_t* Updates_index, + const double breakdown, + double* Slater_inv, + double* determinant) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return qmckl_failwith(context, - QMCKL_NULL_CONTEXT, - "qmckl_sm_naive", - NULL); + return qmckl_failwith( + context, + QMCKL_NULL_CONTEXT, + "qmckl_sm_naive", + NULL); } #ifdef HAVE_HPC @@ -562,26 +581,28 @@ qmckl_exit_code qmckl_sm_naive(const qmckl_context context, } } else { // Updating smaller sub-matrix - return qmckl_sm_naive_hpc(context, - LDS, - Dim, - N_updates, - Updates, - Updates_index, - breakdown, - Slater_inv, - determinant); + return qmckl_sm_naive_hpc( + context, + LDS, + Dim, + N_updates, + Updates, + Updates_index, + breakdown, + Slater_inv, + determinant); } #else - return qmckl_sm_naive_doc(context, - LDS, - Dim, - N_updates, - Updates, - Updates_index, - breakdown, - Slater_inv, - determinant); + return qmckl_sm_naive_doc( + context, + LDS, + Dim, + N_updates, + Updates, + Updates_index, + breakdown, + Slater_inv, + determinant); #endif return QMCKL_FAILURE; @@ -621,6 +642,32 @@ interface end interface #+end_src +#+CALL: generate_f_interface(table=qmckl_sm_naive_args,rettyp=get_value("FRetType"),fname=get_value("Name")) + +#+RESULTS: +#+begin_src f90 :tangle (eval fh_func) :comments org :exports none +interface + integer(c_int32_t) function qmckl_sm_naive & + (context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) & + bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + + integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: LDS + integer (c_int64_t) , intent(in) , value :: Dim + integer (c_int64_t) , intent(in) , value :: N_updates + real (c_double ) , intent(in) :: Updates(N_updates*LDS) + integer (c_int64_t) , intent(in) :: Updates_index(N_updates) + real (c_double ) , intent(in) , value :: breakdown + real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS) + real (c_double ) , intent(inout) :: determinant + + end function qmckl_sm_naive +end interface +#+end_src + #+CALL: generate_f_interface(table=qmckl_sm_naive_args,rettyp=get_value("FRetType"),fname="qmckl_sm_naive_doc") #+RESULTS: @@ -647,32 +694,6 @@ interface end interface #+end_src -#+CALL: generate_f_interface(table=qmckl_sm_naive_args,rettyp=get_value("FRetType"),fname="qmckl_sm_naive") - -#+RESULTS: -#+begin_src f90 :tangle (eval fh_func) :comments org :exports none -interface - integer(c_int32_t) function qmckl_sm_naive & - (context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) & - bind(C) - use, intrinsic :: iso_c_binding - import - implicit none - - integer (c_int64_t) , intent(in) , value :: context - integer (c_int64_t) , intent(in) , value :: LDS - integer (c_int64_t) , intent(in) , value :: Dim - integer (c_int64_t) , intent(in) , value :: N_updates - real (c_double ) , intent(in) :: Updates(N_updates*LDS) - integer (c_int64_t) , intent(in) :: Updates_index(N_updates) - real (c_double ) , intent(in) , value :: breakdown - real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS) - real (c_double ) , intent(inout) :: determinant - - end function qmckl_sm_naive -end interface -#+end_src - *** Performance This function performs best when there is only 1 rank-1 update in the update cycle. It is not useful to use Sherman-Morrison with update splitting for these cycles since splitting @@ -803,27 +824,32 @@ integer function qmckl_sm_splitting_core_doc_f( & updates_index, & breakdown, & s_inv, & - later_updates, & - later_index, & - later, & + later_upds, & + Later_index, & + Later, & determinant) result(info) use qmckl implicit none - integer*8 , intent(in) :: context - integer*8 , intent(in) :: lds, dim - integer*8 , intent(in) :: nupdates - integer*8 , intent(in) :: updates_index(nupdates) - real*8 , intent(in) :: upds(nupdates * lds) - real*8 , intent(in) :: breakdown - real*8 , intent(inout) :: s_inv(dim * lds) - real*8 , intent(inout) :: later_updates(nupdates * lds) - integer*8 , intent(inout) :: later_index(nupdates) - integer*8 , intent(inout) :: later - real*8 , intent(inout) :: determinant + integer*8 , intent(in) :: context + integer*8 , intent(in) :: lds, dim + integer*8 , intent(in) :: nupdates + integer*8 , intent(in) :: updates_index(nupdates) + real*8 , intent(in) :: upds(nupdates * lds) + real*8 , intent(in) :: breakdown + real*8 , intent(inout) :: s_inv(dim * lds) + real*8 , intent(inout) :: determinant + integer*8 , intent(inout) :: Later + integer*8 , intent(inout) :: Later_index(nupdates) + real*8 , intent(inout) :: later_upds(nupdates * lds) - real*8 , dimension(lds, nupdates) :: Updates + real*8 , dimension(nupdates, lds) :: Updates + real*8 , dimension(nupdates, lds) :: Later_updates real*8 , dimension(dim, lds) :: Inverse + real*8 , dimension(dim) :: C + real*8 , dimension(lds) :: D + real*8 :: denominator, idenominator, update + integer*8 :: i, j, l, row info = QMCKL_FAILURE @@ -831,15 +857,65 @@ integer function qmckl_sm_splitting_core_doc_f( & info = QMCKL_INVALID_CONTEXT return endif - + ! Convert 'upds' and 's_inv' into the more easily readable Fortran ! matrices 'Updates' and 'Inverse'. call convert(upds, s_inv, Updates, Inverse, nupdates, lds, dim) - ! YET TO BE IMPLEMENTED - ! Copy updated inverse back to s_inv - call copy_back(Inverse, s_inv, lds, dim) + l = 1; + ! For each update do... + do while (l < nupdates + 1) + + ! Compute C = S^{-1}U(l) + do i = 1, dim + C(i) = 0 + do j = 1, dim + C(i) = C(i) + Inverse(i, j) * Updates(j, l) + end do + end do + + ! Compute denominator = 1 + V(l)^TC + row = updates_index(l) + denominator = 1 + C(row) + + ! If denominator is too close to zero: + ! - Split update in 2 before storing in Later_updates + ! - Split previously computed vector C in 2 + ! - Recompute the denominator + if (abs(denominator) < breakdown) then + do i = 1, dim + Later_updates(i, l) = Updates(i, l) / 2 + C(i) = C(i) / 2 + end do + Later_index(Later) = updates_index(l) + Later = Later + 1 + denominator = 1 + C(row) + end if + + idenominator = 1 / denominator + + ! Update det(S) + determinant = determinant * denominator + + ! selecting column: v_l^T * S_inv + D = Inverse(row, :) + + ! A^{-1} = A^{-1} - C x D / denominator + do i = 1, dim + do j = 1, dim + update = C(i) * D(j) * idenominator + Inverse(i, j) = Inverse(i, j) - update + end do + end do + + l = l + 1 + end do + + ! Copy updated inverse and later updates + ! back to s_inv and later_upds + call copy_back_inv(Inverse, s_inv, lds, dim) + call copy_back_lu(Later_Updates, later_upds, lds, nupdates) info = QMCKL_SUCCESS @@ -855,57 +931,6 @@ for C users and in the module file 'qmckl_f.F90' for Fortran users. #+CALL: generate_c_interface(table=qmckl_sm_splitting_core_args,rettyp=get_value("CRetType"),fname="qmckl_sm_splitting_core_doc") -#+RESULTS: -#+begin_src f90 :tangle (eval f) :comments org :exports none -integer(c_int32_t) function qmckl_sm_splitting_core_doc & - (context, & - LDS, & - Dim, & - N_updates, & - Updates, & - Updates_index, & - breakdown, & - Slater_inv, & - later_updates, & - later_index, & - later, & - determinant) & - bind(C) result(info) - - use, intrinsic :: iso_c_binding - implicit none - - integer (c_int64_t) , intent(in) , value :: context - integer (c_int64_t) , intent(in) , value :: LDS - integer (c_int64_t) , intent(in) , value :: Dim - integer (c_int64_t) , intent(in) , value :: N_updates - real (c_double ) , intent(in) :: Updates(N_updates*LDS) - integer (c_int64_t) , intent(in) :: Updates_index(N_updates) - real (c_double ) , intent(in) , value :: breakdown - real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS) - real (c_double ) , intent(inout) :: later_updates(N_updates*LDS) - integer (c_int64_t) , intent(inout) :: later_index(N_updates) - integer (c_int64_t) , intent(inout) :: later - real (c_double ) , intent(inout) :: determinant - - integer(c_int32_t), external :: qmckl_sm_splitting_core_doc_f - info = qmckl_sm_splitting_core_doc_f & - (context, & - LDS, & - Dim, & - N_updates, & - Updates, & - Updates_index, & - breakdown, & - Slater_inv, & - later_updates, & - later_index, & - later, & - determinant) - -end function qmckl_sm_splitting_core_doc -#+end_src - *** C headers (exposed in qmckl.h) #+CALL: generate_c_header(table=qmckl_sm_splitting_core_args,rettyp=get_value("CRetType"),fname=get_value("Name")) @@ -926,10 +951,29 @@ qmckl_exit_code qmckl_sm_splitting_core ( double* determinant ); #+end_src +#+CALL: generate_c_header(table=qmckl_sm_splitting_core_args,rettyp=get_value("CRetType"),fname="qmckl_sm_splitting_core_hpc") + +#+RESULTS: +#+begin_src c :tangle (eval h_private_func) :comments org +qmckl_exit_code qmckl_sm_splitting_core_hpc ( + const qmckl_context context, + const uint64_t LDS, + const uint64_t Dim, + const uint64_t N_updates, + const double* Updates, + const uint64_t* Updates_index, + const double breakdown, + double* Slater_inv, + double* later_updates, + uint64_t* later_index, + uint64_t* later, + double* determinant ); +#+end_src + #+CALL: generate_c_header(table=qmckl_sm_splitting_core_args,rettyp=get_value("CRetType"),fname="qmckl_sm_splitting_core_doc") #+RESULTS: -#+begin_src c :tangle (eval h_func) :comments org +#+begin_src c :tangle (eval h_func) :no-expand comments org qmckl_exit_code qmckl_sm_splitting_core_doc ( const qmckl_context context, const uint64_t LDS, @@ -1386,7 +1430,8 @@ able to do numerically correct computations, it does not do it in the most effic not be used in real workloads. #+begin_src f90 :tangle (eval f) -integer function qmckl_sm_splitting_doc_f(context, & +integer recursive function qmckl_sm_splitting_doc_f( & + context, & lds, dim, & nupdates, & upds, & @@ -1397,28 +1442,54 @@ integer function qmckl_sm_splitting_doc_f(context, & use qmckl implicit none - integer*8 , intent(in) :: context - integer*8 , intent(in) :: lds, dim - integer*8 , intent(in) :: nupdates - integer*8 , intent(in) :: updates_index(nupdates) - real*8 , intent(in) :: upds(nupdates * lds) - real*8 , intent(in) :: breakdown - real*8 , intent(inout) :: s_inv(dim * lds) - real*8 , intent(inout) :: determinant + integer*8 , intent(in) :: context + integer*8 , intent(in) :: lds, dim + integer*8 , intent(in) :: nupdates + integer*8 , intent(in) :: updates_index(nupdates) + real*8 , intent(in) :: upds(nupdates * lds) + real*8 , intent(in) :: breakdown + real*8 , intent(inout) :: s_inv(dim * lds) + real*8 , intent(inout) :: determinant - real*8 , dimension(lds, nupdates) :: Updates - real*8 , dimension(dim, lds) :: Inverse + integer*8 :: Later + integer*8 , dimension(nupdates) :: Later_index + real*8 , dimension(nupdates * lds) :: Later_updates info = QMCKL_FAILURE - ! Convert 'upds' and 's_inv' into the more easily readable Fortran - ! matrices 'Updates' and 'Inverse'. - call convert(upds, s_inv, Updates, Inverse, nupdates, lds, dim) + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif - ! YET TO BE IMPLEMENTED + Later = 0 + Later_index = 0 + Later_updates = 0 - ! Copy updated inverse back to s_inv - call copy_back(Inverse, s_inv, lds, dim) + info = qmckl_sm_splitting_core_doc_f( & + context, & + lds, dim, & + nupdates, & + upds, & + updates_index, & + breakdown, & + s_inv, & + Later_updates, & + Later_index, & + Later, & + determinant) + + if (Later > 0) then + info = qmckl_sm_splitting_doc_f( & + context, & + lds, dim, & + Later, & + Later_updates, & + Later_index, & + breakdown, & + s_inv, & + determinant) + end if info = QMCKL_SUCCESS From 7a97aa4a77c913a90226e0dc446ab0294929616c Mon Sep 17 00:00:00 2001 From: Francois Coppens Date: Thu, 23 Feb 2023 15:57:35 +0100 Subject: [PATCH 4/8] Fixed Fortran function call bug. --- org/qmckl_sherman_morrison_woodbury.org | 28 +++++++++++++------------ 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/org/qmckl_sherman_morrison_woodbury.org b/org/qmckl_sherman_morrison_woodbury.org index d323cab..c2d1127 100644 --- a/org/qmckl_sherman_morrison_woodbury.org +++ b/org/qmckl_sherman_morrison_woodbury.org @@ -29,7 +29,7 @@ int main() { This is the range that determines the how many high performance kernel instantces will be generated, using the C-function templates defined in the sections below. If the name of the C-function template is called ~qmckl_kernel_{Dim}~, then ~range(K, L+1)~ will results in kernel instances from ~qmckl_kernel_K~ to ~qmckl_kernel_L~. #+NAME:kernel_generator_range #+begin_src python :noweb yes :exports none -range(2, 3) +range(2, 22) #+end_src @@ -244,7 +244,7 @@ integer function qmckl_sm_naive_doc_f(context, & end do ! Copy updated inverse back to s_inv - call copy_back(Inverse, s_inv, lds, dim) + call copy_back_inv(Inverse, s_inv, lds, dim) info = QMCKL_SUCCESS @@ -1451,6 +1451,8 @@ integer recursive function qmckl_sm_splitting_doc_f( & real*8 , intent(inout) :: s_inv(dim * lds) real*8 , intent(inout) :: determinant + integer , external :: qmckl_sm_splitting_core_doc_f + integer*8 :: Later integer*8 , dimension(nupdates) :: Later_index real*8 , dimension(nupdates * lds) :: Later_updates @@ -1690,17 +1692,7 @@ qmckl_exit_code qmckl_sm_splitting( Slater_inv, determinant); #else - // return qmckl_sm_splitting_doc( - // context, - // LDS, - // Dim, - // N_updates, - // Updates, - // Updates_index, - // breakdown, - // Slater_inv, - // determinant); - return qmckl_sm_splitting_hpc( + return qmckl_sm_splitting_doc( context, LDS, Dim, @@ -1710,6 +1702,16 @@ qmckl_exit_code qmckl_sm_splitting( breakdown, Slater_inv, determinant); + // return qmckl_sm_splitting_hpc( + // context, + // LDS, + // Dim, + // N_updates, + // Updates, + // Updates_index, + // breakdown, + // Slater_inv, + // determinant); #endif return QMCKL_SUCCESS; From 8e2674a3b2f9f8062555ed8228848c360b7889bd Mon Sep 17 00:00:00 2001 From: Francois Coppens Date: Thu, 23 Feb 2023 17:16:55 +0100 Subject: [PATCH 5/8] reorder --- org/qmckl_sherman_morrison_woodbury.org | 52 ++++++++++++------------- 1 file changed, 26 insertions(+), 26 deletions(-) diff --git a/org/qmckl_sherman_morrison_woodbury.org b/org/qmckl_sherman_morrison_woodbury.org index c2d1127..831abda 100644 --- a/org/qmckl_sherman_morrison_woodbury.org +++ b/org/qmckl_sherman_morrison_woodbury.org @@ -616,32 +616,6 @@ qmckl_exit_code qmckl_sm_naive(const qmckl_context context, :FRetType: qmckl_exit_code :END: -#+CALL: generate_f_interface(table=qmckl_sm_naive_args,rettyp=get_value("FRetType"),fname="qmckl_sm_naive_hpc") - -#+RESULTS: -#+begin_src f90 :tangle (eval fh_func) :comments org :exports none -interface - integer(c_int32_t) function qmckl_sm_naive_hpc & - (context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) & - bind(C) - use, intrinsic :: iso_c_binding - import - implicit none - - integer (c_int64_t) , intent(in) , value :: context - integer (c_int64_t) , intent(in) , value :: LDS - integer (c_int64_t) , intent(in) , value :: Dim - integer (c_int64_t) , intent(in) , value :: N_updates - real (c_double ) , intent(in) :: Updates(N_updates*LDS) - integer (c_int64_t) , intent(in) :: Updates_index(N_updates) - real (c_double ) , intent(in) , value :: breakdown - real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS) - real (c_double ) , intent(inout) :: determinant - - end function qmckl_sm_naive_hpc -end interface -#+end_src - #+CALL: generate_f_interface(table=qmckl_sm_naive_args,rettyp=get_value("FRetType"),fname=get_value("Name")) #+RESULTS: @@ -668,6 +642,32 @@ interface end interface #+end_src +#+CALL: generate_f_interface(table=qmckl_sm_naive_args,rettyp=get_value("FRetType"),fname="qmckl_sm_naive_hpc") + +#+RESULTS: +#+begin_src f90 :tangle (eval fh_func) :comments org :exports none +interface + integer(c_int32_t) function qmckl_sm_naive_hpc & + (context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) & + bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + + integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: LDS + integer (c_int64_t) , intent(in) , value :: Dim + integer (c_int64_t) , intent(in) , value :: N_updates + real (c_double ) , intent(in) :: Updates(N_updates*LDS) + integer (c_int64_t) , intent(in) :: Updates_index(N_updates) + real (c_double ) , intent(in) , value :: breakdown + real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS) + real (c_double ) , intent(inout) :: determinant + + end function qmckl_sm_naive_hpc +end interface +#+end_src + #+CALL: generate_f_interface(table=qmckl_sm_naive_args,rettyp=get_value("FRetType"),fname="qmckl_sm_naive_doc") #+RESULTS: From 656d2681873204923d6b87c2ec61e91464d290de Mon Sep 17 00:00:00 2001 From: Francois Coppens Date: Sun, 26 Feb 2023 12:26:33 +0100 Subject: [PATCH 6/8] qmckl_sm_splitting_doc kernel works. --- org/qmckl_sherman_morrison_woodbury.org | 125 +++++++++++++++++------- 1 file changed, 88 insertions(+), 37 deletions(-) diff --git a/org/qmckl_sherman_morrison_woodbury.org b/org/qmckl_sherman_morrison_woodbury.org index 831abda..36e62bc 100644 --- a/org/qmckl_sherman_morrison_woodbury.org +++ b/org/qmckl_sherman_morrison_woodbury.org @@ -153,7 +153,7 @@ end subroutine copy_back_inv subroutine copy_back_lu(Later_updates, later_upds, lds, nupdates) implicit none integer*8 , intent(in) :: lds, nupdates - real*8 , intent(in) , dimension(nupdates, lds) :: Later_updates + real*8 , intent(in) , dimension(lds, nupdates) :: Later_updates real*8 , intent(out) :: later_upds(nupdates * lds) integer*8 :: i, j @@ -161,7 +161,7 @@ subroutine copy_back_lu(Later_updates, later_upds, lds, nupdates) ! Copy updated inverse back to s_inv do i = 1, nupdates do j = 1, lds - later_upds((i - 1) * lds + j) = Later_updates(i, j) + later_upds((i - 1) * lds + j) = Later_updates(j, i) end do end do end subroutine copy_back_lu @@ -303,7 +303,7 @@ qmckl_exit_code qmckl_sm_naive ( double* determinant ); #+end_src -#+CALL: generate_c_header(table=qmckl_sm_naive_args,rettyp=get_value("CRetType"),fname="qmckl_sm_naive_hpc") +#+CALL: generate_private_c_header(table=qmckl_sm_naive_args,rettyp=get_value("CRetType"),fname="qmckl_sm_naive_hpc") #+RESULTS: #+begin_src c :tangle (eval h_private_func) :comments org @@ -343,6 +343,7 @@ Common includes and macros used by all the Sherman-Morrison-Woodbury kernels. #include "qmckl.h" #include "config.h" #include "assert.h" +#include "stdio.h" // Order important because // __GNUC__ also set in ICC, ICX and CLANG @@ -835,22 +836,24 @@ integer function qmckl_sm_splitting_core_doc_f( & integer*8 , intent(in) :: lds, dim integer*8 , intent(in) :: nupdates integer*8 , intent(in) :: updates_index(nupdates) - real*8 , intent(in) :: upds(nupdates * lds) + real*8 , intent(in) :: upds(lds * nupdates) real*8 , intent(in) :: breakdown real*8 , intent(inout) :: s_inv(dim * lds) real*8 , intent(inout) :: determinant integer*8 , intent(inout) :: Later integer*8 , intent(inout) :: Later_index(nupdates) - real*8 , intent(inout) :: later_upds(nupdates * lds) + real*8 , intent(inout) :: later_upds(lds * nupdates) - real*8 , dimension(nupdates, lds) :: Updates - real*8 , dimension(nupdates, lds) :: Later_updates + real*8 , dimension(lds, nupdates) :: Updates + real*8 , dimension(lds, nupdates) :: Later_updates real*8 , dimension(dim, lds) :: Inverse real*8 , dimension(dim) :: C real*8 , dimension(lds) :: D real*8 :: denominator, idenominator, update integer*8 :: i, j, l, row + write(*,*) "Entering 'qmckl_sm_splittinig_core_doc_f'" + info = QMCKL_FAILURE if (context == QMCKL_NULL_CONTEXT) then @@ -888,7 +891,7 @@ integer function qmckl_sm_splitting_core_doc_f( & Later_updates(i, l) = Updates(i, l) / 2 C(i) = C(i) / 2 end do - Later_index(Later) = updates_index(l) + Later_index(Later + 1) = updates_index(l) Later = Later + 1 denominator = 1 + C(row) end if @@ -919,6 +922,8 @@ integer function qmckl_sm_splitting_core_doc_f( & info = QMCKL_SUCCESS + write(*,*) "Leaving 'qmckl_sm_splittinig_core_doc_f'" + end function qmckl_sm_splitting_core_doc_f #+end_src @@ -931,6 +936,62 @@ for C users and in the module file 'qmckl_f.F90' for Fortran users. #+CALL: generate_c_interface(table=qmckl_sm_splitting_core_args,rettyp=get_value("CRetType"),fname="qmckl_sm_splitting_core_doc") +#+RESULTS: +#+begin_src f90 :tangle (eval f) :comments org :exports none +integer(c_int32_t) function qmckl_sm_splitting_core_doc & + (context, & + LDS, & + Dim, & + N_updates, & + Updates, & + Updates_index, & + breakdown, & + Slater_inv, & + later_updates, & + later_index, & + later, & + determinant) & + bind(C) result(info) + + use, intrinsic :: iso_c_binding + implicit none + + integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: LDS + integer (c_int64_t) , intent(in) , value :: Dim + integer (c_int64_t) , intent(in) , value :: N_updates + real (c_double ) , intent(in) :: Updates(N_updates*LDS) + integer (c_int64_t) , intent(in) :: Updates_index(N_updates) + real (c_double ) , intent(in) , value :: breakdown + real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS) + real (c_double ) , intent(inout) :: later_updates(N_updates*LDS) + integer (c_int64_t) , intent(inout) :: later_index(N_updates) + integer (c_int64_t) , intent(inout) :: later + real (c_double ) , intent(inout) :: determinant + + integer(c_int32_t), external :: qmckl_sm_splitting_core_doc_f + + write(*,*) "Entering 'qmckl_sm_splittinig_core_doc'" + + info = qmckl_sm_splitting_core_doc_f & + (context, & + LDS, & + Dim, & + N_updates, & + Updates, & + Updates_index, & + breakdown, & + Slater_inv, & + later_updates, & + later_index, & + later, & + determinant) + + write(*,*) "Leaving 'qmckl_sm_splittinig_core_doc'" + +end function qmckl_sm_splitting_core_doc +#+end_src + *** C headers (exposed in qmckl.h) #+CALL: generate_c_header(table=qmckl_sm_splitting_core_args,rettyp=get_value("CRetType"),fname=get_value("Name")) @@ -1221,7 +1282,7 @@ qmckl_exit_code qmckl_sm_splitting_core( if (LDS == (1+(Dim-1)/SIMD_LENGTH)*SIMD_LENGTH) { // Most cases switch (Dim) { <> - case default: { + default: { assert(0 == 1 && "TEMPLATE NOT IMPLEMENTED!"); break; } @@ -1243,20 +1304,7 @@ qmckl_exit_code qmckl_sm_splitting_core( determinant); } #else - // return qmckl_sm_splitting_core_doc( - // context, - // LDS, - // Dim, - // N_updates, - // Updates, - // Updates_index, - // breakdown, - // Slater_inv, - // later_updates, - // later_index, - // later, - // determinant); - return qmckl_sm_splitting_core_hpc( + return qmckl_sm_splitting_core_doc( context, LDS, Dim, @@ -1446,7 +1494,7 @@ integer recursive function qmckl_sm_splitting_doc_f( & integer*8 , intent(in) :: lds, dim integer*8 , intent(in) :: nupdates integer*8 , intent(in) :: updates_index(nupdates) - real*8 , intent(in) :: upds(nupdates * lds) + real*8 , intent(in) :: upds(lds * nupdates) real*8 , intent(in) :: breakdown real*8 , intent(inout) :: s_inv(dim * lds) real*8 , intent(inout) :: determinant @@ -1455,7 +1503,9 @@ integer recursive function qmckl_sm_splitting_doc_f( & integer*8 :: Later integer*8 , dimension(nupdates) :: Later_index - real*8 , dimension(nupdates * lds) :: Later_updates + real*8 , dimension(lds * nupdates) :: Later_updates + + write(*,*) "Entering 'qmckl_sm_splitting_doc_f'" info = QMCKL_FAILURE @@ -1495,6 +1545,8 @@ integer recursive function qmckl_sm_splitting_doc_f( & info = QMCKL_SUCCESS + write(*,*) "Leaving 'qmckl_sm_splitting_doc_f'" + end function qmckl_sm_splitting_doc_f #+end_src @@ -1520,16 +1572,21 @@ integer(c_int32_t) function qmckl_sm_splitting_doc & integer (c_int64_t) , intent(in) , value :: LDS integer (c_int64_t) , intent(in) , value :: Dim integer (c_int64_t) , intent(in) , value :: N_updates - real (c_double ) , intent(in) :: Updates(N_updates*LDS) + real (c_double ) , intent(in) :: Updates(LDS*N_updates) integer (c_int64_t) , intent(in) :: Updates_index(N_updates) real (c_double ) , intent(in) , value :: breakdown real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS) real (c_double ) , intent(inout) :: determinant integer(c_int32_t), external :: qmckl_sm_splitting_doc_f + + write(*,*) "Entering 'qmckl_sm_splitting_doc'" + info = qmckl_sm_splitting_doc_f & (context, LDS, Dim, N_updates, Updates, Updates_index, breakdown, Slater_inv, determinant) + write(*,*) "Leaving 'qmckl_sm_splitting_doc'" + end function qmckl_sm_splitting_doc #+end_src @@ -1672,6 +1729,8 @@ qmckl_exit_code qmckl_sm_splitting( const double breakdown, double* Slater_inv, double* determinant) { + + printf("Entering 'qmckl_sm_splitting'\n"); if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith( @@ -1680,7 +1739,7 @@ qmckl_exit_code qmckl_sm_splitting( "qmckl_sm_splitting", NULL); } - #ifdef HAS_HPC + #ifdef HAVE_HPC return qmckl_sm_splitting_hpc( context, LDS, @@ -1702,18 +1761,10 @@ qmckl_exit_code qmckl_sm_splitting( breakdown, Slater_inv, determinant); - // return qmckl_sm_splitting_hpc( - // context, - // LDS, - // Dim, - // N_updates, - // Updates, - // Updates_index, - // breakdown, - // Slater_inv, - // determinant); #endif + printf("Leaving 'qmckl_sm_splitting'\n"); + return QMCKL_SUCCESS; } #+end_src From 1640eb60f97b1574af22dae982ee06edb942fdd9 Mon Sep 17 00:00:00 2001 From: Francois Coppens Date: Mon, 27 Feb 2023 11:29:01 +0100 Subject: [PATCH 7/8] Changed argument order. --- org/qmckl_sherman_morrison_woodbury.org | 23 +++++++++-------------- 1 file changed, 9 insertions(+), 14 deletions(-) diff --git a/org/qmckl_sherman_morrison_woodbury.org b/org/qmckl_sherman_morrison_woodbury.org index 36e62bc..efe10fc 100644 --- a/org/qmckl_sherman_morrison_woodbury.org +++ b/org/qmckl_sherman_morrison_woodbury.org @@ -790,11 +790,11 @@ If the determinant is passed it will only be partially updated if there were any | ~LDS~ | ~uint64_t~ | in | Leading dimension of Slater_inv | | ~Dim~ | ~uint64_t~ | in | Dimension of Slater_inv | | ~N_updates~ | ~uint64_t~ | in | Number of rank-1 updates to be applied to Slater_inv | -| ~Updates~ | ~double[N_updates*LDS]~ | in | Array containing the rank-1 updates | +| ~Updates~ | ~double[LDS*N_updates]~ | in | Array containing the rank-1 updates | | ~Updates_index~ | ~uint64_t[N_updates]~ | in | Array containing positions of the rank-1 updates | | ~breakdown~ | ~double~ | in | Break-down parameter on which to fail or not | | ~Slater_inv~ | ~double[Dim*LDS]~ | inout | Array containing the inverse Slater-matrix | -| ~later_updates~ | ~double[N_updates*LDS]~ | inout | Array containing the split updates for later | +| ~later_updates~ | ~double[LDS*N_updates]~ | inout | Array containing the split updates for later | | ~later_index~ | ~uint64_t[N_updates]~ | inout | Array containing the positions of the split updates for later | | ~later~ | ~uint64_t~ | inout | Number of split updates for later | | ~determinant~ | ~double~ | inout | Determinant of the Slater-matrix | @@ -928,11 +928,11 @@ end function qmckl_sm_splitting_core_doc_f #+end_src **** C interface to the pedagogical kernel (not directly exposed) -The following Fortran function ~qmckl_sm_splitting_core_doc~ makes sure -that the pedagogical kernel ~qmckl_sm_splitting_core_doc_f~, written in -Fortran, can be called from C using the ~ISO_C_BINDING~. The Fortran function -~qmckl_sm_splitting_core_doc~ will be exposed in the header file 'qmckl.h' -for C users and in the module file 'qmckl_f.F90' for Fortran users. +The function ~qmckl_sm_splitting_core_doc~ makes sure that +~qmckl_sm_splitting_core_doc_f~ can be called from C using the +~ISO_C_BINDING~. Function ~qmckl_sm_splitting_core_doc~ will be +exposed in ~qmckl.h~ and ~qmckl_f.F90~, but +~qmckl_sm_splitting_core_doc_f~ will not. #+CALL: generate_c_interface(table=qmckl_sm_splitting_core_args,rettyp=get_value("CRetType"),fname="qmckl_sm_splitting_core_doc") @@ -960,19 +960,16 @@ integer(c_int32_t) function qmckl_sm_splitting_core_doc & integer (c_int64_t) , intent(in) , value :: LDS integer (c_int64_t) , intent(in) , value :: Dim integer (c_int64_t) , intent(in) , value :: N_updates - real (c_double ) , intent(in) :: Updates(N_updates*LDS) + real (c_double ) , intent(in) :: Updates(LDS*N_updates) integer (c_int64_t) , intent(in) :: Updates_index(N_updates) real (c_double ) , intent(in) , value :: breakdown real (c_double ) , intent(inout) :: Slater_inv(Dim*LDS) - real (c_double ) , intent(inout) :: later_updates(N_updates*LDS) + real (c_double ) , intent(inout) :: later_updates(LDS*N_updates) integer (c_int64_t) , intent(inout) :: later_index(N_updates) integer (c_int64_t) , intent(inout) :: later real (c_double ) , intent(inout) :: determinant integer(c_int32_t), external :: qmckl_sm_splitting_core_doc_f - - write(*,*) "Entering 'qmckl_sm_splittinig_core_doc'" - info = qmckl_sm_splitting_core_doc_f & (context, & LDS, & @@ -987,8 +984,6 @@ integer(c_int32_t) function qmckl_sm_splitting_core_doc & later, & determinant) - write(*,*) "Leaving 'qmckl_sm_splittinig_core_doc'" - end function qmckl_sm_splitting_core_doc #+end_src From b01c7c306be0ef4d8c544654c59dc41bc8ae6a3f Mon Sep 17 00:00:00 2001 From: Francois Coppens Date: Mon, 27 Feb 2023 12:04:04 +0100 Subject: [PATCH 8/8] Done with splitting doc. --- org/qmckl_sherman_morrison_woodbury.org | 24 +++++++++++------------- 1 file changed, 11 insertions(+), 13 deletions(-) diff --git a/org/qmckl_sherman_morrison_woodbury.org b/org/qmckl_sherman_morrison_woodbury.org index efe10fc..4ddd949 100644 --- a/org/qmckl_sherman_morrison_woodbury.org +++ b/org/qmckl_sherman_morrison_woodbury.org @@ -371,7 +371,7 @@ include situations where one wants to apply updates to a square submatrix of the full matrix. It takes advantage of memory aligned data and assumes no data dependencies inside the loops. The loops are fully vectorised whenever ~Dim~ is an integer -multiple of ~SIMD_LEGTH~. +multiple of ~SIMD_LENGTH~. #+begin_src c :tangle (eval c) :comments org qmckl_exit_code qmckl_sm_naive_hpc( const qmckl_context context, @@ -525,9 +525,9 @@ text=""" result = [] for Dim in <>: Dim=str(Dim) - result.append(text.replace("{Dim}",Dim) ) + result.append(text.replace("{Dim}",Dim)) -return '\n'.join(result) +return ''.join(result) #+end_src Python script that generated C switch cases that call individual kernel instances. @@ -541,14 +541,13 @@ case {Dim}: Updates_index, breakdown, Slater_inv, - determinant); -""" + determinant);""" result = [] for Dim in <>: Dim=str(Dim) - result.append(text.replace("{Dim}",Dim) ) + result.append(text.replace("{Dim}",Dim)) -return '\n'.join(result) +return ''.join(result) #+end_src #+begin_src c :tangle (eval c) :comments org :noweb yes @@ -1225,7 +1224,7 @@ for Dim in <>: Dim=str(Dim) result.append(text.replace("{Dim}",Dim) ) -return '\n'.join(result) +return ''.join(result) #+end_src #+NAME:slagel_splitting_switch-case_generator @@ -1244,14 +1243,13 @@ case {Dim}: { later, determinant); break; -} -""" +}""" result = [] for Dim in <>: Dim=str(Dim) result.append(text.replace("{Dim}",Dim) ) -return '\n'.join(result) +return ''.join(result) #+end_src #+begin_src c :tangle (eval c) :comments org :noweb yes @@ -1603,10 +1601,10 @@ qmckl_exit_code qmckl_sm_splitting ( double* determinant ); #+end_src -#+CALL: generate_c_header(table=qmckl_sm_splitting_args,rettyp=get_value("CRetType"),fname="qmckl_sm_splitting_hpc") +#+CALL: generate_private_c_header(table=qmckl_sm_splitting_args,rettyp=get_value("CRetType"),fname="qmckl_sm_splitting_hpc") #+RESULTS: -#+begin_src c :tangle (eval h_func) :comments org +#+begin_src c :tangle (eval h_private_func) :comments org qmckl_exit_code qmckl_sm_splitting_hpc ( const qmckl_context context, const uint64_t LDS,