diff --git a/org/qmckl_sherman_morrison_woodbury.org b/org/qmckl_sherman_morrison_woodbury.org index ef21226..bebfc70 100644 --- a/org/qmckl_sherman_morrison_woodbury.org +++ b/org/qmckl_sherman_morrison_woodbury.org @@ -21,8 +21,6 @@ Low- and high-level functions that use the Sherman-Morrison and Woodbury matrix #define THRESHOLD 1e-3 #endif -#include "qmckl.h" - int main() { qmckl_context context; context = qmckl_context_create(); @@ -31,7 +29,7 @@ int main() { * Sherman-Morrison Helper Functions -[TODO: FMJC] Add doc +Helper functions that are used by the Sherman-Morrison-Woodbury kernels. These functions should only be used in the context of these kernels. ** ~qmckl_sherman_morrison_threshold~ :PROPERTIES: @@ -40,7 +38,7 @@ int main() { :FRetType: double precision :END: -[TODO: FMJC] Add doc +This function is used to set the threshold value that is used in the kernels to determine if a matrix is invertable or not. In the Sherman-Morrison kernels this is determined by comparing the denominator in the Sherman-Morrison formula to the value set in threshold. If the value is smaller than the threshold value it means the matrix is not invertable. In the Woodbury kernels the threshold value is compared with the value of the determinant of the update matrix. #+NAME: qmckl_sherman_morrison_threshold_args | double | thresh | out | Threshold | @@ -87,9 +85,9 @@ end function qmckl_sherman_morrison_threshold_f // Sherman-Morrison-Woodbury break-down threshold static double qmckl_shreman_morrison_threshold(double* const threshold) { *threshold = THRESHOLD; -#ifdef DEBUG - std::cerr << "Break-down threshold set to: " << threshold << std::endl; -#endif +// #ifdef DEBUG +// std::cerr << "Break-down threshold set to: " << threshold << std::endl; +// #endif } #+end_src @@ -216,9 +214,9 @@ qmckl_exit_code qmckl_sherman_morrison_c(const qmckl_context context, const double* Updates, const uint64_t* Updates_index, double * Slater_inv) { -#ifdef DEBUG - std::cerr << "Called qmckl_sherman_morrison with " << N_updates << " updates" << std::endl; -#endif +// #ifdef DEBUG +// std::cerr << "Called qmckl_sherman_morrison with " << N_updates << " updates" << std::endl; +// #endif double C[Dim]; double D[Dim]; @@ -239,8 +237,10 @@ qmckl_exit_code qmckl_sherman_morrison_c(const qmckl_context context, // Denominator double den = 1 + C[Updates_index[l] - 1]; - if (fabs(den) < threshold) { - return false; + double const thresh; + qmckl_sherman_morrison_threshold(&thresh); + if (fabs(den) < thresh) { + return QMCKL_FAILURE; } double iden = 1 / den;