Move verify gradient functions to another helper file

vgvassilev · Nov 20, 2024 · fa5fdde · fa5fdde
1 parent 4bc4b79
commit fa5fdde
Show file tree

Hide file tree

Showing 2 changed files with 124 additions and 105 deletions.
diff --git a/demos/CUDA/BlackScholes/BlackScholes.cu b/demos/CUDA/BlackScholes/BlackScholes.cu
@@ -35,10 +35,12 @@
  * DISCLAIMER: The following file has been slightly modified to ensure
  * compatibility with Clad and to serve as a Clad demo. Specifically, parts of
  * the original `main` function have been moved to a separate function to use
- * `clad::gradient` on. Furthermore, original print statements have been removed
- * and new helper functions are now included in the file to verify the
- * gradient's results. The original file is available in NVIDIA's cuda-samples
- * repository on GitHub.
+ * `clad::gradient` on. Furthermore, Clad cannot clone printf statements so some
+ * original print statements have been ommitted. The same applies to the
+ * checkCudaErrors function.
+ * New helper functions are included in another file and invoked here to verify
+ * the gradient's results. The original file is available in NVIDIA's
+ * cuda-samples repository on GitHub.
  *
  * Relevant documentation regarding the problem at hand can be found in NVIDIA's
  * cuda-samples repository. Using Clad, we compute some of the Greeks
@@ -53,6 +55,7 @@
 
 #include <helper_cuda.h> // helper functions CUDA error checking and initialization
 #include <helper_functions.h> // helper functions for string parsing
+#include <helper_grad_verify.h>
 
 ////////////////////////////////////////////////////////////////////////////////
 // Process an array of optN options on CPU
@@ -77,17 +80,18 @@ float RandFloat(float low, float high) {
   return (1.0f - t) * low + t * high;
 }
 
+// This section is included in the helper_grad_verify.h file
 ////////////////////////////////////////////////////////////////////////////////
 // Data configuration
 ////////////////////////////////////////////////////////////////////////////////
-const int OPT_N = 4000000;
-const int NUM_ITERATIONS = 512;
+// const int OPT_N = 4000000;
+// const int NUM_ITERATIONS = 512;
 
-const int OPT_SZ = OPT_N * sizeof(float);
-const float RISKFREE = 0.02f;
-const float VOLATILITY = 0.30f;
+// const int OPT_SZ = OPT_N * sizeof(float);
+// const float RISKFREE = 0.02f;
+// const float VOLATILITY = 0.30f;
 
-#define DIV_UP(a, b) (((a) + (b) - 1) / (b))
+// #define DIV_UP(a, b) (((a) + (b) - 1) / (b))
 
 ////////////////////////////////////////////////////////////////////////////////
 // Main program
@@ -139,101 +143,6 @@ void launch(float* h_CallResultCPU, float* h_CallResultGPU,
   cudaFree(d_CallResult);
 }
 
-double d1(double S, double X, double T) {
-  return (log(S / X) + (RISKFREE + 0.5 * VOLATILITY * VOLATILITY) * T) /
-         (VOLATILITY * sqrt(T));
-}
-
-double N_prime(double d) {
-  const double RSQRT2PI =
-      0.39894228040143267793994605993438; // 1 / sqrt(2 * PI)
-  return RSQRT2PI * exp(-0.5 * d * d);
-}
-
-enum Greek { Delta, dX, Theta };
-
-enum OptionType { Call, Put };
-
-template <OptionType opt> const char* getNameofOpt() {
-  if constexpr (opt == Call)
-    return "Call";
-  if constexpr (opt == Put)
-    return "Put";
-}
-
-template <Greek greek> const char* getNameOfGreek() {
-  if constexpr (greek == Delta)
-    return "Delta";
-  if constexpr (greek == dX)
-    return "dStrike";
-  if constexpr (greek == Theta)
-    return "Theta";
-}
-
-template <OptionType opt, Greek greek>
-void computeL1norm(float* S, float* X, float* T, float* d) {
-  double delta, ref, sum_delta, sum_ref;
-  sum_delta = 0;
-  sum_ref = 0;
-  for (int i = 0; i < OPT_N; i++) {
-    if constexpr (opt == Call) {
-      if constexpr (greek == Delta) {
-        double d1_val = d1(S[i], X[i], T[i]);
-        ref = CND(d1_val);
-      } else if constexpr (greek == dX) {
-        double T_val = T[i];
-        double d1_val = d1(S[i], X[i], T_val);
-        double d2_val = d1_val - VOLATILITY * sqrt(T_val);
-        double expRT = exp(-RISKFREE * T_val);
-        ref = -expRT * CND(d2_val);
-      } else if constexpr (greek == Theta) {
-        double S_val = S[i], X_val = X[i], T_val = T[i];
-        double d1_val = d1(S_val, X_val, T_val);
-        double d2_val = d1_val - VOLATILITY * sqrt(T_val);
-        double expRT = exp(-RISKFREE * T_val);
-        ref = (S_val * N_prime(d1_val) * VOLATILITY) / (2 * sqrt(T_val)) +
-              RISKFREE * X_val * expRT *
-                  CND(d2_val); // theta is with respect to t, so -theta is the
-                               // approximation of the derivative with respect
-                               // to T
-      }
-    } else if constexpr (opt == Put) {
-      if constexpr (greek == Delta) {
-        double d1_val = d1(S[i], X[i], T[i]);
-        ref = CND(d1_val) - 1.0;
-      } else if constexpr (greek == dX) {
-        double T_val = T[i];
-        double d1_val = d1(S[i], X[i], T_val);
-        double d2_val = d1_val - VOLATILITY * sqrt(T_val);
-        double expRT = exp(-RISKFREE * T_val);
-        ref = expRT * CND(-d2_val);
-      } else if constexpr (greek == Theta) {
-        double S_val = S[i], X_val = X[i], T_val = T[i];
-        double d1_val = d1(S_val, X_val, T_val);
-        double d2_val = d1_val - VOLATILITY * sqrt(T_val);
-        double expRT = exp(-RISKFREE * T_val);
-        ref = (S_val * N_prime(d1_val) * VOLATILITY) / (2 * sqrt(T_val)) -
-              RISKFREE * X_val * expRT * CND(-d2_val);
-      }
-    }
-
-    delta = fabs(d[i] - ref);
-    sum_delta += delta;
-    sum_ref += fabs(ref);
-  }
-
-  double L1norm = sum_delta / sum_ref;
-  printf("L1norm of %s for %s option = %E\n", getNameOfGreek<greek>(),
-         getNameofOpt<opt>(), L1norm);
-  if (L1norm > 1e-5) {
-    printf(
-        "Gradient test failed: Difference between %s's computed and "
-        "approximated theoretical values for %s option is larger than expected",
-        getNameOfGreek<greek>(), getNameofOpt<opt>());
-    exit(EXIT_FAILURE);
-  }
-}
-
 int main(int argc, char** argv) {
   // Start logs
   printf("[%s] - Starting...\n", argv[0]);

diff --git a/demos/CUDA/BlackScholes/helper/helper_grad_verify.h b/demos/CUDA/BlackScholes/helper/helper_grad_verify.h
@@ -0,0 +1,110 @@
+#include <cmath>
+
+extern "C" double CND(double d);
+
+////////////////////////////////////////////////////////////////////////////////
+// Data configuration
+////////////////////////////////////////////////////////////////////////////////
+const int OPT_N = 4000000;
+const int NUM_ITERATIONS = 512;
+
+const int OPT_SZ = OPT_N * sizeof(float);
+const float RISKFREE = 0.02f;
+const float VOLATILITY = 0.30f;
+
+#define DIV_UP(a, b) (((a) + (b)-1) / (b))
+
+double d1(double S, double X, double T) {
+  return (log(S / X) + (RISKFREE + 0.5 * VOLATILITY * VOLATILITY) * T) /
+         (VOLATILITY * sqrt(T));
+}
+
+double N_prime(double d) {
+  const double RSQRT2PI =
+      0.39894228040143267793994605993438; // 1 / sqrt(2 * PI)
+  return RSQRT2PI * exp(-0.5 * d * d);
+}
+
+enum Greek { Delta, dX, Theta };
+
+enum OptionType { Call, Put };
+
+template <OptionType opt> const char* getNameofOpt() {
+  if constexpr (opt == Call)
+    return "Call";
+  if constexpr (opt == Put)
+    return "Put";
+}
+
+template <Greek greek> const char* getNameOfGreek() {
+  if constexpr (greek == Delta)
+    return "Delta";
+  if constexpr (greek == dX)
+    return "dStrike";
+  if constexpr (greek == Theta)
+    return "Theta";
+}
+
+template <OptionType opt, Greek greek>
+void computeL1norm(float* S, float* X, float* T, float* d) {
+  double delta, ref, sum_delta, sum_ref;
+  sum_delta = 0;
+  sum_ref = 0;
+  for (int i = 0; i < OPT_N; i++) {
+    if constexpr (opt == Call) {
+      if constexpr (greek == Delta) {
+        double d1_val = d1(S[i], X[i], T[i]);
+        ref = CND(d1_val);
+      } else if constexpr (greek == dX) {
+        double T_val = T[i];
+        double d1_val = d1(S[i], X[i], T_val);
+        double d2_val = d1_val - VOLATILITY * sqrt(T_val);
+        double expRT = exp(-RISKFREE * T_val);
+        ref = -expRT * CND(d2_val);
+      } else if constexpr (greek == Theta) {
+        double S_val = S[i], X_val = X[i], T_val = T[i];
+        double d1_val = d1(S_val, X_val, T_val);
+        double d2_val = d1_val - VOLATILITY * sqrt(T_val);
+        double expRT = exp(-RISKFREE * T_val);
+        ref = (S_val * N_prime(d1_val) * VOLATILITY) / (2 * sqrt(T_val)) +
+              RISKFREE * X_val * expRT *
+                  CND(d2_val); // theta is with respect to t, so -theta is the
+                               // approximation of the derivative with respect
+                               // to T
+      }
+    } else if constexpr (opt == Put) {
+      if constexpr (greek == Delta) {
+        double d1_val = d1(S[i], X[i], T[i]);
+        ref = CND(d1_val) - 1.0;
+      } else if constexpr (greek == dX) {
+        double T_val = T[i];
+        double d1_val = d1(S[i], X[i], T_val);
+        double d2_val = d1_val - VOLATILITY * sqrt(T_val);
+        double expRT = exp(-RISKFREE * T_val);
+        ref = expRT * CND(-d2_val);
+      } else if constexpr (greek == Theta) {
+        double S_val = S[i], X_val = X[i], T_val = T[i];
+        double d1_val = d1(S_val, X_val, T_val);
+        double d2_val = d1_val - VOLATILITY * sqrt(T_val);
+        double expRT = exp(-RISKFREE * T_val);
+        ref = (S_val * N_prime(d1_val) * VOLATILITY) / (2 * sqrt(T_val)) -
+              RISKFREE * X_val * expRT * CND(-d2_val);
+      }
+    }
+
+    delta = fabs(d[i] - ref);
+    sum_delta += delta;
+    sum_ref += fabs(ref);
+  }
+
+  double L1norm = sum_delta / sum_ref;
+  printf("L1norm of %s for %s option = %E\n", getNameOfGreek<greek>(),
+         getNameofOpt<opt>(), L1norm);
+  if (L1norm > 1e-5) {
+    printf(
+        "Gradient test failed: Difference between %s's computed and "
+        "approximated theoretical values for %s option is larger than expected",
+        getNameOfGreek<greek>(), getNameofOpt<opt>());
+    exit(EXIT_FAILURE);
+  }
+}