gauss legendre integration

build.bat

@@ -22,6 +22,9 @@ set libiompdll=%libiompdll_path%\%libiompdll_name%
 
 set Sources=../src/main.c
 
+
+
+IF NOT EXIST .\out mkdir .\out
 IF NOT EXIST .\build mkdir .\build
 pushd .\build
 

src/main.c

@@ -1,4 +1,5 @@
 #include <stdlib.h>
+#include <math.h>
 
 #define ENABLE_LOGGING 1
 #if ENABLE_LOGGING
@@ -70,9 +71,66 @@ struct Atom
 };
 
 
+typedef struct GaussLegendre GaussLegendre;
+struct GaussLegendre
+{
+  U32 order;
+  F64 *weights;
+  F64 *abscissae;
+};
+
+
 //~ Globals
 global Grid g_grid = {0};
 global BSplineCtx g_bspline_ctx = {0};
+global GaussLegendre g_gauss_legendre = {0};
+
+
+function void
+set_up_gauss_legendre_points(Arena* arena)
+{
+  U32 GL_order = 7;
+
+  if(GL_order == 7)
+  {
+    F64 quadrature_weights[7] = {
+            0.41795918367346938,
+            0.38183005050511894,
+            0.38183005050511894,
+            0.27970539148927666,
+            0.27970539148927666,
+            0.12948496616886969,
+            0.12948496616886969
+    };
+
+    F64 legendre_abscissae[7] = {
+            0,
+            0.405845151377397,
+            -0.40584515137739,
+            -0.74153118559939,
+            0.741531185599394,
+            -0.94910791234275,
+            0.949107912342758
+    }; 
+
+    g_gauss_legendre.order = GL_order;
+    g_gauss_legendre.weights = PushArray(arena, F64, GL_order);
+    g_gauss_legendre.abscissae = PushArray(arena, F64, GL_order);
+    for(U32 i = 0; i < GL_order; i++) 
+    {
+      g_gauss_legendre.weights[i] = quadrature_weights[i];
+      g_gauss_legendre.abscissae[i] = legendre_abscissae[i];
+    }
+
+  }
+
+}
+
+
+
+
+
+
 
 function void
 write_array_binary_F64(String8 path_to_file, F64 *values, U32 array_size) 
@@ -340,7 +398,8 @@ function void bspline_things()
 
 /* Auxiliary routine: printing a matrix */
 function void
-print_matrix_Z64( char* desc, int m, int n, Z64* a, int lda ) {
+print_matrix_Z64( char* desc, int m, int n, Z64* a, int lda ) 
+{
   ArenaTemp scratch = scratch_get(0,0);
   int i, j;
   String8 newline = str8_lit("\n");
@@ -358,6 +417,77 @@ print_matrix_Z64( char* desc, int m, int n, Z64* a, int lda ) {
   }
 }
 
+typedef F64 (*func_F64)(F64);
+
+function F64
+cos2(F64 x) 
+{
+  return cos(x)*cos(x);
+}
+
+function F64
+gq_integration(F64 a, F64 b, func_F64 func) 
+{
+
+  F64 prefac = 0.5*(b-a);
+  F64 gq_sum = 0.0;
+
+  for(U32 i = 0; i < g_gauss_legendre.order; i++) 
+  {
+    F64 w = g_gauss_legendre.weights[i];
+    F64 z = g_gauss_legendre.abscissae[i];
+    F64 shift = (z*prefac)+((a+b)*0.5);
+    F64 funct_value_at_shift = func(shift);
+    gq_sum = gq_sum + prefac*w*funct_value_at_shift;
+  }
+
+  return gq_sum;
+
+}
+
+function void
+test_gauss_legendre()
+{
+
+  // Test GL qudrature by integrating cos^2 from 0 to 2pi, it should equal pi.
+  {
+    ArenaTemp scratch = scratch_get(0,0);
+
+    F64 pi = 4.0*atan(1.0);
+    String8 pi_out = str8_pushf(scratch.arena, "Pi from standard library 4.0*atan(1.0) = %.16f \n", pi);
+    LOG(pi_out.str);
+
+    F64 a = 0.0;
+    F64 b = 2.0*pi;
+    F64 gq_sum_single_interval = gq_integration(a, b, cos2);
+
+    String8 out = str8_pushf(scratch.arena, "Integration result for a single interval: %.16f, Reference: %.16f \n", gq_sum_single_interval, pi); 
+    LOG(out.str);
+
+    // Test with several smaller intervals instead.
+    F64 aa[10] = {0.0};
+    F64 bb[10] = {0.0};
+    F64 delta = b/10.0;
+    for(U32 i = 0; i < 10; i++) 
+    {
+      aa[i] = i*delta;
+      bb[i] = (i+1)*delta;
+      String8 intervals = str8_pushf(scratch.arena, "%i, a=%f, b=%f \n", i, aa[i], bb[i]);
+      LOG(intervals.str);
+    }
+
+    F64 gq_sum_several_intervals = 0.0;
+    for(U32 i = 0; i < 10; i++) 
+    {
+      gq_sum_several_intervals += gq_integration(aa[i], bb[i], cos2);
+
+    }
+
+    out = str8_pushf(scratch.arena, "Integration result for ten intervals: %.16f, Reference: %.16f \n", gq_sum_several_intervals, pi); 
+    LOG(out.str);
+  }
+}
+
 
 
 function void EntryPoint(void) 
@@ -367,7 +497,28 @@ function void EntryPoint(void)
   OS_InitGfxReceipt os_gfx_receipt = OS_gfx_init(os_receipt);
 
   Arena *arena = m_make_arena();
- 
+
+  set_up_gauss_legendre_points(arena);
+
+  
+
+
+
+
+
+}
+
+
+
+
+function void mkl_things(void) 
+{
+
+  OS_InitReceipt os_receipt = OS_init();
+  OS_InitGfxReceipt os_gfx_receipt = OS_gfx_init(os_receipt);
+
+  Arena *arena = m_make_arena();
+
   U32 N = 4;
   Z64 *main_A = PushArray(arena, Z64, N*N);
   main_A[0]  = (Z64){-3.84,  2.25};