Added threshold function. Tests still dont compile. #25

org/qmckl_sherman_morrison_woodbury.org

@@ -28,6 +28,120 @@ int main() {
   context = qmckl_context_create();
 
   #+end_src
+
+* Sherman-Morrison Helper Functions
+
+[TODO: FMJC] Add doc
+
+** ~qmckl_sherman_morrison_threshold~
+   :PROPERTIES:
+   :Name:     qmckl_sherman_morrison_threshold
+   :CRetType: double
+   :FRetType: double precision
+   :END:
+
+[TODO: FMJC] Add doc
+
+   #+NAME: qmckl_sherman_morrison_threshold_args
+   | double | thresh | out | Threshold |
+
+*** Requirements
+
+      Add description of the input variables. (see for e.g. qmckl_distance.org)
+
+*** C header
+
+    #+CALL: generate_c_header(table=qmckl_sherman_morrison_threshold_args,rettyp=get_value("CRetType"),fname=get_value("Name"))
+
+    #+RESULTS:
+    #+begin_src c :tangle (eval h_func) :comments org
+// Sherman-Morrison-Woodbury break-down threshold
+#ifndef THRESHOLD
+#define THRESHOLD 1e-3
+#endif
+
+    qmckl_exit_code qmckl_sherman_morrison_threshold (
+	  double* const thresh ); 
+    #+end_src
+
+
+*** Source Fortran
+
+    #+begin_src f90 :tangle (eval f)
+integer function qmckl_sherman_morrison_threshold_f(thresh) result(info) 
+  use qmckl
+  implicit none
+  real*8     , intent(inout)  :: thresh
+  !logical, external :: qmckl_sherman_morrison_f
+  info = qmckl_sherman_morrison_threshold(thresh)
+end function qmckl_sherman_morrison_threshold_f
+    #+end_src
+
+*** Source C
+
+    #+begin_src c :tangle (eval c) :comments org
+#include <stdbool.h>
+#include <math.h>
+#include "qmckl.h"
+
+// 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
+}
+    #+end_src
+
+
+*** Performance
+
+** C interface                                                     :noexport:
+
+   #+CALL: generate_c_interface(table=qmckl_sherman_morrison_threshold_args,rettyp=get_value("FRetType"),fname=get_value("Name"))
+
+   #+RESULTS:
+   #+begin_src f90 :tangle (eval f) :comments org :exports none
+   integer(c_int32_t) function qmckl_sherman_morrison_threshold &
+       (thresh) &
+       bind(C) result(info)
+
+     use, intrinsic :: iso_c_binding
+     implicit none
+
+     real    (c_double ) , intent(out)         :: thresh
+
+     integer(c_int32_t), external :: qmckl_sherman_morrison_threshold_f
+     info = qmckl_sherman_morrison_threshold_f &
+	    (thresh)
+
+   end function qmckl_sherman_morrison_threshold
+   #+end_src
+
+   #+CALL: generate_f_interface(table=qmckl_sherman_morrison_threshold_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_sherman_morrison_threshold &
+	 (thresh) &
+	 bind(C)
+       use, intrinsic :: iso_c_binding
+       import
+       implicit none
+
+       real    (c_double ) , intent(out)         :: thresh
+
+     end function qmckl_sherman_morrison_threshold
+   end interface
+   #+end_src
+
+*** Test                                                           :noexport:
+
+[TODO: FMJC] Write tests for the Sherman-Morrison part.
+
+
+
   
 * Naïve   Sherman-Morrison
   
@@ -64,12 +178,6 @@ int main() {
 
     #+RESULTS:
     #+begin_src c :tangle (eval h_func) :comments org
-// Sherman-Morrison-Woodbury break-down threshold
-#ifndef THRESHOLD
-#define THRESHOLD 1e-3
-#endif
-static double threshold();
-
     qmckl_exit_code qmckl_sherman_morrison_c (
 	  const qmckl_context context,
 	  const uint64_t Dim,
@@ -102,15 +210,6 @@ end function qmckl_sherman_morrison_f
 #include <stdbool.h>
 #include "qmckl.h"
 
-// Sherman-Morrison-Woodbury break-down threshold
-static double threshold() {
-  const double threshold = THRESHOLD;
-#ifdef DEBUG
-  std::cerr << "Break-down threshold set to: " << threshold << std::endl;
-#endif
-  return threshold;
-}
-
 qmckl_exit_code qmckl_sherman_morrison_c(const qmckl_context context, 
                                 const uint64_t Dim,
                                 const uint64_t N_updates,
@@ -124,6 +223,9 @@ qmckl_exit_code qmckl_sherman_morrison_c(const qmckl_context context,
   double C[Dim];
   double D[Dim];
 
+  double threshold = 0.0;
+  qmckl_sherman_morrison_threshold(&threshold);
+
   unsigned int l = 0;
   // For each update
   while (l < N_updates) {
@@ -137,7 +239,7 @@ qmckl_exit_code qmckl_sherman_morrison_c(const qmckl_context context,
 
     // Denominator
     double den = 1 + C[Updates_index[l] - 1];
-    if (fabs(den) < threshold()) {
+    if (fabs(den) < threshold) {
       return false;
     }
     double iden = 1 / den;